- Email: [email protected]
Numerical analysis of phase behavior during rapid decompression of rich natural gases
Accepted Manuscript Title: Numerical analysis of phase behavior during rapid decompression of rich natural gases Author: Evgeniy Burlutskiy PII: DOI: Reference:
S0957-5820(13)00014-1 http://dx.doi.org/doi:10.1016/j.psep.2013.03.001 PSEP 351
To appear in:
Process Safety and Environment Protection
Received date: Revised date: Accepted date:
21-11-2012 5-3-2013 12-3-2013
Please cite this article as: Burlutskiy, E., Numerical analysis of phase behavior during rapid decompression of rich natural gases, Process Safety and Environment Protection (2013), http://dx.doi.org/10.1016/j.psep.2013.03.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical analysis of phase behavior during rapid decompression of rich natural gases
Evgeniy Burlutskiy
ip t
A*STAR Institute of High Performance Computing (IHPC), 1 Fusionopolis Way, 138632 Singapore
cr
Tel: (65) 6419 1386 Fax: (65) 6467 4350 e-mail: [email protected]
us
Highlights 1D transient model of gas-liquid two-phase flow in a shock tube is developed
an
Approach to model the condensation in a compressible gas mixture is proposed The model is validated on the experimental data on RGD in rich natural gases
M
Physical reasons form the “condensation plateau” is investigated numerically
d
Abstract
te
The effect of the condensation process on the gas and liquid phase behavior during rapid decompression of rich natural gases is studied in the paper numerically. A one-dimensional
Ac ce p
mathematical model of transient thermal two-phase flow of compressible multi-component natural gas mixture and liquid phase in a shock tube is developed. The set of mass, momentum and enthalpy conservation equations is solved for the gas and liquid phases. The approach to model a liquid condensation process during rapid decompression of rich natural gas mixture is proposed. The mass transfer between the gas and the liquid is taken into account by introducing the appropriate terms into the governing equations. Thermo-physical properties of multicomponent natural gas mixture are calculated by solving the Equation of State (EOS) in the form of the Soave-Redlich-Kwong (SRK-EOS) model. The proposed liquid condensation model is integrated into the GDP model. A simple case of GDP model, where the liquid was not considered, was extensively validated on base and dry natural gases. The proposed two-phase
1
Page 1 of 38
model is validated against the experiments where the decompression wave speed was measured in rich natural gases at low temperature. It shows a good agreement with the experimental data.
KEY WORDS: Mathematical modeling; multiphase flow; Mass transfer; Condensation; Gases;
cr
ip t
Rapid decompression.
us
Introduction
an
Gas processing plants as well as basic fossil-fuel customers move further away from source production fields and carbon capture and storage (CCS) places. Rich and lean natural
M
gases and carbon dioxide (CO2) mixtures are transported by pipeline systems which are very long in length. Fractures in pipes caused by ruptures bring much more severe consequences
d
nowadays because of high operating pressure. A crack may easy propagate and destroy the
te
entire pipeline. The fracture propagation control in pipeline systems is usually based on the employment of the Battelle two-curve method (BTCM) developed by the Battelle Columbus
Ac ce p
Laboratories in order to determine the minimum of fracture arrest toughness (Eiber et. al., 1993, 2004). The crack propagation is arrested if the decompression wave speed in the gas is quicker than the fracture propagation velocity in the pipeline wall. Analytical models such as GASDECOM (Eiber et. al., 1993) are usually used by gas pipeline engineers in order to calculate the thermo-physical properties of natural gases. The program predicts the decompression wave speed in natural gases with a reasonably high level of accuracy (Botros et. al., 2003, 2007) using different equation of state (EOS) models. The effect of gas-to-wall friction is not taken into account because the fluid transport equations are not considered. The GASDECOM's predictive performance is poor if a relatively narrow shock tube (where the inner diameter is less than 5-10 cm) having a significant roughness of the internal surface is considered. The GASDECOM model does not work in the case of gas-liquid two-phase flows 2
Page 2 of 38
and if liquid condensation which is usually observed in the case of rich natural gases and carbon dioxide (CO2) mixtures is occurred. The effect of liquid on the decompression process may be accounted for in GASDECOM manually by introducing an average phase properties but the accuracy of such solutions is poor. However, more detailed information about time history
ip t
values and spatial distribution of basic flow parameters is also important. A mathematical model of fast transient multi-component fluid flows in a shock tube is highly desirable in this case in
cr
order to simulate a complex fluid dynamics. The problem became multiphase if the
us
condensation is occurred in the case of rich natural gases and carbon dioxide (CO2) mixtures.
an
Experimental investigation of single- and two-phase fluid dynamics in a shock tube is much easy nowadays because the accuracy of experimental instruments and the computer power
M
are significantly increased. First full scale experiment on the decompression (Tam and Cowley, 1988) were performed by British Petroleum and Shell Oil on the Isle of Grain (UK). A
d
horizontal pipe having a length of 100 meters and a diameter of 15 cm was used. The initial
te
pressure was varied from 0.8 MPa to 2.1 MPa in this study. The LPG mixture consisting of 95% mole propane and 5% mole butane was used as the fluid phase. Wide experimental
Ac ce p
measurements of the decompression wave speed in rich, base and dry natural gas mixtures were conducted by Trans Canada Pipe Lines (Botros et. al., 2003, 2004, 2007, 2010a, 2010b). Most of the measurements were made by using small-diameter shock tubes, where the friction influences on the flow behavior much stronger as compared to large-diameter pipes. The effect of shock tube diameter, gas mixture composition and fluid thermo-physical properties was examined experimentally. Pressure values were varied in those studies in the range between 10 MPa and 37 MPa (Botros et. al., 2003, 2010a). The temperature was varied from normal values to low (Botros et. al., 2003, 2007, 2010a). Time history data collections of basic flow parameters were obtained experimentally for rich natural gases in Botros et. al. (2007).
3
Page 3 of 38
Mathematical modeling of multi-phase fluid dynamics in a shock tube having natural gases at high pressure was also conducted. A simple method of solving the problem of multiphase blowdown flow from short pipes was proposed by Chen et. al. (1995a, 1995b). The set of two-phase flow governing equations was solved for multi-component fluid mixtures. The
ip t
proposed model was validated on the experiments of rapid decompression of simple gases. Numerical predictions showed a reasonable fit to the experiments of Tam and Cowley (1988). A
cr
numerical model of fluid dynamics following full-bore rupture in pipe systems having multicomponent fluid mixtures was proposed in Mahgerefteh et. al. (2006a, 2006b). The set of
us
governing equations was solved by using the method of characteristics (Zucrow and Hoffman,
an
1976). The model assumes a homogeneous equilibrium of the fluid flow, where the constituent phases are at thermal and mechanical equilibrium. The model focuses on simulation of complex
M
large-scale pipeline systems having bends, branches and couplings. The model was validated through the experiments of Tam and Cowley (1988) and a hypothetical pipeline system having a
d
length of 25 km successfully. A condensation process was not considered in those studies. A
te
one-dimensional code OLGA developed by SPT-group (Benediksen et. al., 1991) was examined in (Botros et. al., 2007) as a numerical solver for the problem of rapid decompression of dry and
Ac ce p
base natural gas mixtures. The attempt was un-successful and simulations showed a significant over-prediction of the pressure time history and a poor comparison (Burlutskiy 2013b) with the experimental data.
More general approach was applied in Burlutskiy (2012a, 2013a) for single gas and gas-
liquid flows in a shock tube (GDP model). The set of governing equations was solved for the gas and liquid phases together with taking into account of the Equation of State (EOS) model. The pressure time evolution in a shock tube for different axial elevations was predicted directly using the mass balance approach (Burlutskiy, 2013a, 2013b). The fracture propagation control was also successfully performed using the proposed GDP model (Burlutskiy, 2012b) in the frame of the Battelle two-curve method (BTCM), where two curves such as the crack velocity and the 4
Page 4 of 38
decompression wave speed were numerically simulated and compared with each other. Gasliquid system where the liquid is represented by the film phase and the condensation process were not considered in those works but the GDP model showed a great potential in modelling of
ip t
the liquid condensation process.
Numerical investigation of the condensation effect on the gas and liquid phase’s
cr
behavior during rapid decompression of rich natural gases is presented in the paper. The analysis is performed using a one-dimensional mathematical model of transient thermal two-phase flow
us
of compressible multi-component gas mixture and liquid film phase in a shock tube. The
an
approach to model the process of liquid condensation during rapid decompression of rich natural gases as well as the numerical algorithm of solving the problem is proposed in the paper. The
M
proposed model is validated against the experiments (Botros et. al., 2007) where the
te
d
decompression wave speed in rich natural gases was measured.
Ac ce p
One-dimensional mathematical model of gas-liquid two-phase flow in a shock tube
A one-dimensional transient mathematical model of compressible thermal gas-liquid
two-phase flow in a shock tube is developed. The set of mass, momentum and enthalpy conservation equations is solved for the gas and liquid phases. The liquid phase is assumed to be in the form of liquid film because the internal pipe diameter in the considered case is very small (49.325 mm) as compared to industrial pipeline diameters (about 1 m). The mass transfer has a form of fast condensation process from the gaseous phase into the liquid film phase. The natural gas mixture is considered as the gas core phase in this system. The set of governing equations in general form is written as (Wallis, 1969):
5
Page 5 of 38
G G G G U G G F t z
(1)
F F F F U F t z
(2)
G F
2 G G U G G G U G P G t z z
F
P z
R GF
G F U G F
R G F
R G Wall
R F Wall
(3)
cr
F F U F F F U 2F t z
G F U G F
ip t
an
P (U UG )2 P q G F G F G F F UF z 2 t R G F ( U F U G ) R F Wall U F
(5)
(6)
M
F F h F F F U Fh F t z
us
P G G h G G G U G h G P (U U G )2 q G F G F G F G UG t z z 2 t R G F ( U F U G ) R G Wall U G
(4)
Here, G , F are the volume fraction of gas mixture and liquid film, respectively; G ,
d
F are the density of the gas and the film; UG , UF are the velocity of the gas and the film,
te
respectively; P is the total pressure; R G Wall , R F Wall are friction terms, which represent the
Ac ce p
gas-to-wall and the film-to-wall friction, respectively; R G F is the term, which represents the friction between the gas and film phases; h G , h F are the enthalpy of the gas and the film; q G F is the term, which represents a heat exchange between the gas phase and the liquid film phase; G F is the term, which represents a condensation from the gaseous phase into the liquid phase;
t is the time; z is the axial co-ordinate. Last terms in the enthalpy conservation equations (5)-(6) represent additional contribution into the system’s enthalpy due to the friction force. It is assumed that the velocity of gas fraction, which leaves the gaseous phase during the condensation process and passes into the liquid state, is equal to the velocity of the liquid fraction (i.e. U G F U F ) during the condensation.
6
Page 6 of 38
The mass transfer term which describes a liquid condensation process is represented by the following correlation:
i MC S
(7)
ip t
G F
Here, M C is the condensation rate. The mass transfer term (7) is a universal correlation which is
cr
dependent from the shock tube inner diameter (i.e. i D pipe ; S D 2pipe / 4 ); i is the
us
perimeter of the gas core phase at core-annular flow regime (i.e. the gas core-to-liquid film interface); S is the cross-sectional area of the pipe. Therefore, it can be easy extrapolated into
an
other cases. The condensation rate ( M C ) is a constant in the proposed study and its value will be determined numerically through a better fit to the experimental data. A dependence of the
M
condensation rate ( M C ) from the various flow parameters such as the gas mixture composition, temperature, pressure and others will be investigated in other studies. The friction term, which
G Wall S
(8)
Ac ce p
R G Wall
te
d
represents gas-to-wall friction, is written in the Blasius form (Blasius, 1913):
G Wall
2 , Re G 1600 G Wall 64 / Re G G Wall G U G , 0.25 8 G Wall 0.316 / Re G , Re G 1600
ReG G UG DpipeG / G
(9)
(10)
Here, is the perimeter of the pipe; G Wall is the gas-to-wall friction term; G Wall is
the friction coefficient, D pipe is the diameter of the pipe; G is the gas viscosity. The term, which describes the friction between the film and the wall, is written in the Taitel and Dukler form (Taitel and Dukler, 1976):
7
Page 7 of 38
R F Wall
F Wall S
F Wall F Wall F U 2F
(11) , Re F 1600 F Wall 64 / Re F F Wall 0.046 / Re 0F.2 , Re F 1600
,
(12)
ReF F UFDpipeF / F
ip t
(13)
cr
Here, F Wall is the liquid film-to-wall friction term; F Wall is the friction coefficient,
F is the viscosity of the liquid. The friction between the gas and the film is also written in the
us
Blasius form (Blasius, 1913):
i G F S
G F
, Re G 1600 G F 64 / Re G G F G U G U F ( U G U F ) , 0 . 25 8 G F 0.316 / Re G , Re G 1600
(14)
(15)
d
M
an
R GF
te
Here, G F is the gas-to-film friction term; G F is the friction coefficient. The heat
Ac ce p
transfer between the gas phase and the liquid film phase is written as (Wallis, 1969):
qGF
i i G Nu G F (T TF ) S 4 S(1 F )
(16)
0.5 0.33 is the Here, CP G is the heat capacity of the gas phase; NuG F 2 0.6 ReG F PrG
Nusselt number of the flow; G is the coefficient of heat conductivity of the gas phase.
Thermo-physical properties of natural gas mixtures
8
Page 8 of 38
Thermo-physical fluid properties are modeled by solving of the Equations of State (EOS) in the form of the Soave-Redlich-Kwong (SRK-EOS) model (Soave, 1979). It has been chosen due to better performance as compared to all other EOS models. The set of equations and
a V V b
N N
a zi z j a i a j 1 kij i 1j1
(17)
;
N
b z i bi i 1
2
an
R 2 TC2 i T 1 m i 1 a i 0.42748 PC i TC i
RTC i PC i
M
bi 0.08664
cr
V b
us
RT
P
ip t
correlations (SRK-EOS) may be written as (Soave, 1979):
(19)
(20)
(21)
d
m i 0 . 48 1 . 574 i 0 . 176 i2
(18)
te
Here, V is the volume of the gas mixture; N is the number of components in the gas
Ac ce p
mixture; T is the temperature of the gas mixture; R is the universal gas constant; i is the acentric factor of the component i; PC i , TC i are critical values of the pressure and temperature, correspondently; zi is the mole fraction of the component i. The compressibility factor (Z) of gas mixture is calculated from the following equation (Soave, 1979):
Z3 Z2 A B B2 Z AB 0 A
aP R 2T 2
P G
;
B
(22)
bP RT
(23)
R TZ M
(24)
9
Page 9 of 38
Here, M is the molar mass of a substance. The viscosity of gas mixture is calculated by using the Lee-Gonzales-Eakin (LGE) correlation (Lee et. al., 1966).
ip t
Solving algorithm of two-phase problem
cr
The algorithm of solving of the set of one-dimensional transient governing equations of the fluid mixture flow in a pipe is based on the Tri-Diagonal Matrix Algorithm (TDMA), also
us
known as the Thomas algorithm (Patankar, 1980). It is a simplified form of Gaussian
an
elimination that can be used to solve tri-diagonal systems of equations. The computational pipe is divided into N parts, which are called the control volumes. The set of unsteady governing
M
equations is transformed into the standard form of the discrete analog of the tri-diagonal system (Patankar, 1980) by using the fully implicit numerical scheme and it is calculated at each control
d
volume of the computational pipe. The equation is reduces to the steady state discretization
te
equation if the time step goes to infinity.
Ac ce p
The gas and liquid initial velocities are set to zero ( U G i = 0 m/s; U L i = 0 m/s; i=1, ..., N;
N is the number of control volumes of the pipe) in each control volume of the computational pipe before the decompression. The pressure and the temperature of the gas mixture are set to specified initial values. The pressure outside the tube is set to the atmospheric ( PN Pin ; PN 1 = 0.1 MPa; PN >> PN 1 ; N+1 is a ghost cell). The initial volume fraction values of the liquid are set to zero. The initial density of the gas mixture is calculated by using the SRK-EOS model. This is the direct way of using the SRK-EOS model, where the initial density of the gas mixture is a function of the initial pressure and temperature (fig. 1). The following assumptions are made in order to model the compressible nature of the process:
10
Page 10 of 38
The gas mixture density is assumed to be constant when the mass and momentum conservation equations are solved for both phases. A change in the gas density is determined after when the set of governing equations is solved for the considered time step.
The time step of
2 10 5
ip t
The liquid phase density is a constant during the decompression process in the study. sec is selected for simulations in order to provide a continuous
us
cr
change of the gas density.
Main iteration loop starts from the momentum transfer predictions (fig. 1). Velocity
an
values of the gas and liquid film phases are calculated by solving the momentum conservation equations of each phase. The mass continuity equations are solved for both phases after solving
M
the momentum transfer. The momentum and mass continuity equations are solved for the liquid film phase if the liquid is already occurred. The volume fraction of the gas is conceived as a
d
local variable for one time step. It is set to 1 (one) at each control volume of the pipe before
te
solving the momentum and mass continuity equations. When the time step is calculated, a new density of the gas mixture is determined. The old gas density values from the previous time step
Ac ce p
are multiplied on the gas volume fraction values from new time step in order to account for the momentum change for the considered time step. Further, the gas volume fraction values are set to 1 (one) at each control volume for a new time step calculations again. The condensation process is taken into account by introducing the appropriate terms into the governing equations. The mass balance between the gas and the liquid phase for the selected control volume is considered, keeping the total mass of both phases in the selected control volume constant for the considered time step.
When all changes in the gas density is taken into account, new pressure values are determined thereafter. These values are calculated by using the SRK-EOS model determining the inverse correspondence, which may be written schematically as (fig. 1): 11
Page 11 of 38
new old Pnew f (G ,T )
(25)
So, new pressure values are determined at each control volume of the pipe. New
ip t
temperature values for both phases are calculated thereafter by solving the enthalpy conservation equations for the gas and liquid phases. Thus, the density, pressure, temperature and velocity
cr
values are new. This procedure is repeated again for a new time step if the total calculation time
us
is smaller than selected prediction time. The mass transfer is started if the pressure and temperature values of the gas mixture are lower than the critical condensation values (i.e.
an
Pi PCOND and Ti TCOND ) which are dependent from gas mixture composition. The condition
of liquid occurring in the shock tube is checked at every time step before solving the set of
M
governing equations. The proposed model was developed under the research project “Fast Transient Flows of Multi-Component Natural Gas Mixture and Liquid in Pipelines” in
d
PETROSOFT-DC. It was implemented into the FORTRAN computer code named the “Gas
Ac ce p
te
Decompression Program” (GDP code).
Numerical simulation of rapid decompression of rich natural gas mixture
A simple case of the proposed GDP model, where the condensation was not considered,
has been validated against the experimental data on rapid decompression of a shock tube having dry and base natural gases (Botros et. al., 2003) and presented in Burlutskiy (2013a, 2013b). Single-phase flows of natural gases at thermal equilibrium were considered. The proposed twophase model, where the condensation process is taken into account and integrated into the GDP code, is validated against the decompression wave speed measurements were made by Trans Canada Pipe Lines (TCPL) at TCPL Gas Dynamic Test Facility in Didsbury, Alberta, Canada (Botros et. al., 2007) using the same experimental facility as in the case of dry and base natural 12
Page 12 of 38
gases (Burlutskiy, 2013a, 2013b). Fig. 2 shows a photo of the real experimental shock tube facility. The main test section is a shock tube, which has a total length of 172 m. The internal pipe diameter is 49.325 mm. It has a roughness better than 1 micro-meter (i.e. less than 1.0 10 6 6
m). Considering the roughness-to-diameter ratio ( / Dpipe 20.3 10 ;
1.0 10 6
m;
ip t
Dpipe 49.325103 m), the internal surface of the shock tube is determined as the “smooth” tube
cr
surface (Botros et. al., 2010b). The classical Blasius correlation (8)-(10) for the smooth pipes is allowed in this case. A rupture disc was placed at one end of the experimental shock tube, which
us
was upon rupturing. A decompression wave propagates up into the pressurized test section. High frequency responses Pressure Transducers (PT) were mounted into the tube part, which was
an
most close to the rupture end of the shock tube, in order to capture the time history of the expansion fan (Botros et. al., 2007). Distances between the pressure transducers and the rupture
M
disc are shown in table 1. The schematic of the decompression tube is shown in Fig. 3. Four
te
wave speed in those experiments.
d
pressure-time tracers (PT-2, PT-3, PT-4 and PT-6) were used to determinate the decompression
Ac ce p
Fig. 2 and Fig. 3 Table 1, 2
The computational decompression pipe has a length of 15 m. The set of governing
equations is solved on the mesh having 1000 grid nodes on the axial direction. The inner pipe diameter is 49.325 mm. One end of the shock tube is selected to be the closed end. The rupture disc condition is modeled in the other end of the pipe. The initial pressure ( Pin ), temperature ( Tin ), and natural gas mixture composition are shown in table 2 (Botros et. al., 2007). The considered rich natural gas mixture is in the gaseous state if the pressure is in the range between 8 and 10 MPa, and the temperature is between -5 o C and -10 o C . The gas mixture is expanded towards the open end of the shock tube if the rupture disc is opened. The pressure and the 13
Page 13 of 38
temperature of the gas are dropped down due to gas expansion. The mixture which is in the gaseous state at the beginning goes to two-phase region of the pressure-temperature envelope if the pressure is less than 8 MPa and the temperature is less than -10 o C . Below those critical P and T values, the mixture is transformed from single phase system (gas) to two-phase system
ip t
(gas-liquid). Those critical pressure and temperature parameters are dependent from the natural gas mixture composition. Thermo-physical properties of the considered two-phase gas-liquid
cr
system are calculated by employing the SRK-EOS model together with the flash algorithm
us
(Meulbroek 2002) in the conditions of two-phase region of pressure-temperature envelope. The liquid phase density is determined numerically using the SRK-EOS model and the flash kg / m 3 .
an
algorithm and it is equal to 350
M
Fig. 4
d
A collection of measured (Botros et. al., 2007) and predicted temperature time history
te
values at PT-2 and PT-8 locations are shown in fig. 4. The temperature was not collected at other considered PT’s (PT-8 and PT-10) at the experimental measurements. Experimental points
Ac ce p
are shown as symbols here. Broken curves represent calculated values. Predicted temperature time history values are in a good agreement with the measurements in most part of the pressure range. A decrease rate of predicted time-temperature curves is stronger as compared to measured points for pipe elevations near the rupture end (PT-2). The slope of the predicted curve at PT-8 pipe location is much better correlate with the experimental data. Hence, predicted decompressed temperature values (i.e. the level after the decompression) are in a good correlation with measured temperature minimum for every of the considered PT locations. The proposed model simulates the temperature field of the gas mixture with high order of agreement with measured values.
Fig. 5 14
Page 14 of 38
A comparison between measured (Botros et. al., 2007) and predicted pressure time history values at PT-8 and PT-10 locations is shown in fig. 5. Experimental points are shown as symbols. A continues line and broken curves represent calculated values. Measured pressure
ip t
values are decreased at each PT location of the shock tube if the decompression in the gas mixture is started at those PT locations. A decrease of the pressure is more rapid if the
cr
considered PT location closer to the rupture end of the pipe, and is slower for the PT’s far away from the rupture end due to strong effect of gas-to-wall friction. Hence, predicted pressure time
us
history values are in a good agreement with the measurements in most part of the pressure range.
an
The slope of the calculated time-pressure curves is in a good correlation with the experimental values. Moreover, predicted pressure values at the considered PT locations reach measured
M
pressure minimum after the decompression successfully.
te
d
Fig. 6
A correlation between the decompression wave speed and the pressure has to be
Ac ce p
determined in order to perform the Battelle two-curve analysis. The decompression wave speed is calculated by using predicted pressure time history values at PT-2, PT-3, PT-4 and PT-6 locations. Fig. 6 shows measured (Botros et. al., 2007) and predicted decompression wave speed values as a function of the pressure. Pressure values are normalized on the initial pressure in the shock tube before the decompression ( P / Pin ). The initial pressure is equal to 9.949 MPa in this study. Symbols correspond to measured values (Botros et. al., 2007). Curves are numerical simulations performed by using the proposed model, where the liquid condensation is taken into account. The condensation rate is determined numerically through a better fit to the experiments where the decompression wave speed in rich natural gas mixture was measured. The mass transfer effects on the decompression wave speed when the condensation process is occurred (i.e. the pressure ratio between 0.7 and 0.8 in fig. 6). Different values of the condensation rate 15
Page 15 of 38
( M C ) are examined numerically. Fig. 6 shows predictions of the decompression wave speed if the mass transfer rate ( M C ) is equal to 6.125, 12.25 and 49 (kg/m2/s). The best fit to measured values is observed (Fig. 6) if the condensation rate value is equal to 12.25 (kg/m2/s). Real values of the condensation rate were not measured experimentally. Therefore, the condensation rate
ip t
( M C ) equal to 12.25 (kg/m2/s) is chosen to be used in this study because it provides the best fit to measurements. A rapid decrease in the decompression wave speed is not predicted by the
cr
model if the mass transfer between phases is not taken into account or the condensation rate is
an
Fig. 7
us
too small.
M
Fig. 7 shows a development of the liquid phase front propagation up into the pressurized shock tube section predicted by the proposed model. It is compared with predicted lines which
d
represent the propagation of gas decompression wave front for few specified pressure values.
te
The pipe axial elevation in this figure shows the distance from the rupture end of the shock tube. The decompression wave in rich natural gas mixture propagates almost linearly if the pressure is
Ac ce p
above the critical condensation values (i.e. 9.8 MPa; 9.0 MPa; 8.0 MPa lines in fig. 7). The liquid film propagation shows a non-linear behavior in the beginning. The propagation velocity of the liquid phase front is significantly decreased during first 5-6 milliseconds as compared to the beginning. The liquid phase effects on the development of gas decompression wave propagation if pressure values below the critical (i.e. 7.0 MPa; 6.0 MPa lines in fig. 7) during first 2-3 milliseconds due to phases re-distribution and interaction between gas and liquid phases. Such non-linearity in the decompression wave propagation in rich natural gas mixture cannot be predicted directly by using of analytical models.
Fig. 8
16
Page 16 of 38
A flat plateau is observed in both pressure time history curves (fig. 5). This “condensation plateau” is observed if the condensation process in natural gases is occurred only and it is investigated in more details below. For that reason a few more figures which contain time history values at three different specified tube elevations are considered. Distances between
ip t
those elevations and the rupture disc are 0.1 m, 0.5 m and 1 m. Predicted time history values of a slip between velocities of the liquid film and the gas phase at those elevations are shown in fig. 8
cr
and analyzed. The slip is increased at every of those locations because the gas is accelerated when it is released. It consists from the gas only at the beginning because the liquid is occurred
us
later. The slip is increased until a certain level (15-20 m/s) and has a local maximum here. The
an
slip is decreased thereafter because the liquid is condensed and accelerated by the gas, and has a local minimum at the end of the condensation process (i.e. t = 2 ms for z = 0.1 m; t = 4 ms for = 0.5 m; t = 8 ms for z = 1 m). The condensation front keeps propagating up into the
M
z
pressurized pipe section. The decompression continues at those locations due to gas release.
d
However, the gas mixture loses much more kinetic energy than before by accelerating the liquid
te
phase. For the most remote locations ( z = 1 m) gas-to-liquid friction effects much stronger
Ac ce p
(fig. 8) because of a bigger distance from the rupture end, which is occupied by the liquid.
Fig. 9 and Fig. 10
Fig. 9 shows predicted velocity values of the gas and liquid phases which are collected at
the same pipe elevations (0.1 m, 0.5 m and 1 m). Fig. 10 shows pressure time evolution of the gas mixture at those elevations also. Both figures show that the “condensation plateau” in the pressure time history (fig 10, 7 MPa < P < 8 MPa) is bigger for the most distant axial elevation ( z = 1 m). The velocity of the gas mixture is accelerated from 0 m/s to 20 m/s for the elevation of 1 m due to gas expansion at the beginning (fig. 9,
z
= 1 m, 3 ms < t < 5 ms). A pressure
decrease for the considered pipe elevation of 1 m is slower as compared to the elevation of 0.1 m 17
Page 17 of 38
because the friction effects much stronger due to longer distance from the rupture end. A condensation is not occurred yet during this time period at the elevation of 1 m. The liquid condensation front is propagated up from the elevation of 0.6 m to the elevation of 0.8 m during this time period (fig. 7). The velocity of the gas does not show the acceleration and is almost
ip t
constant during the time between 4.5 and 7 ms. Gas expansion in this case is extremely slow due to strong increase of the distance where the liquid condensation is occurred already, and gas-to-
cr
liquid friction contribution, respectively. For the pipe elevation of 0.5 m the effect of the “condensation plateau” is weaker because the distance from the rupture end is shorter. For the
us
pipe elevation of 0.1 m the effect of the “condensation plateau” is absent because the friction is
an
very small in this case. Therefore, high friction between the gas and the liquid phase (gas-toliquid friction) forms the “condensation plateau” in the pressure. The condensation process is
M
started thereafter at the elevation of 1 m when the pressure and the temperature are dropped down and the critical pressure and temperature values are reached (fig. 9,
z
= 1 m, liquid, t = 7
z
= 1 m, liquid, 7 ms < t < 8.5 ms). A pressure decrease is weaker
te
mechanisms act both (fig. 9,
d
ms). The pressure is decreased again because the condensation and momentum transfer
again starting from the time of 8.5 ms because the only momentum transfer (gas release) acts on
Ac ce p
the gas decompression here. Fig. 9 and 10 show together that an increase in gas-to-liquid friction downstream of the considered pipe elevation, due to liquid front propagation toward the pressurized pipe section, forms the “condensation plateau” in the pressure time history.
The decompression in rich natural gas mixtures below the critical condensation
parameters is determined by the action of two mechanisms such as the gas release due to the pressure drop (the momentum transfer) and the liquid condensation process (the mass transfer). A rapid decrease in the decompression wave speed at condensation is dependent from the “condensation plateau” in the pressure time history, which is used to determined those values. Numerical analysis shows that the “condensation plateau” in the pressure time history values is formed due to the following reasons: 18
Page 18 of 38
First, the gas mixture leaves the considered pipe elevation and it is accelerated itself due to gas release. The liquid phase is occurred in the downstream due to the condensation and propagates up into the pressurized pipe section starting from the rupture end. The gas
ip t
interacts with the liquid in the downstream by losing the kinetic energy because it accelerates the liquid phase in the downstream. Such decrease of the gas kinetic energy
cr
in the downstream due to increase of gas-to-liquid friction (the momentum transfer) slows down of the gas acceleration by keeping the gas velocity constant. It reduces the
an
“plateau” in the pressure time history values.
us
pressure drop between the considered pipe elevation and the downstream, and forms the
M
At the later stage, the pressure and the temperature are dropped down enough and the critical pressure and temperature values are reached. The liquid phase is condensed in the
d
considered elevation and the decompression in this location is determined by the action
te
of two mechanisms such as the gas release (momentum transfer) and the condensation process (mass transfer) at the same time together. Therefore, a pressure decrease is
Ac ce p
significantly quicker at the considered elevation as compared to usual decrease rate caused by the momentum transfer only. The pressure decrease is slow again when the condensation is finished.
The effect of the liquid phase on the gas mixture flow behavior is much stronger if small-
diameter pipes and low temperatures are considered. The volume fraction of the liquid phase is much bigger in this case and the friction force effects on the flow behavior much stronger.
Conclusions
19
Page 19 of 38
The effect of the condensation process on the gas and liquid phase behavior during rapid decompression of rich natural gases is analyzed numerically. A one-dimensional mathematical model of transient thermal two-phase flows of compressible multi-component gas mixture and liquid phase in a shock tube is developed. The set of mass, momentum and enthalpy
ip t
conservation equations is solved for the gas and liquid phases. The approach to model the process of liquid condensation during rapid decompression of a compressible rich natural gas
cr
mixture as well as the numerical algorithm of solving the problem is proposed in the paper. The mass transfer process between the gas and the liquid is taken into account by introducing the
us
appropriate terms into the set of governing equations. The proposed model is validated on the
an
experiments where the decompression wave speed was measured in a shock tube having rich
M
natural gases at low temperature and it shows a good agreement with the experimental data.
Numerical analysis shows a non-linear development of the liquid phase front
d
propagation. The propagation speed of liquid phase front is significantly decreased during first
te
few milliseconds as compared to the propagation start. The liquid phase influences on the development of gas decompression wave propagation if the pressure values below the critical
Ac ce p
condensation values during first few milliseconds due to phases re-distribution and frictional interaction between the gas and liquid phases. Such non-linearity in the decompression wave development in natural gas mixtures cannot be predicted directly by using of analytical models.
Numerical analysis shows that the “condensation plateau” in the decompression wave
speed is taken place due to action of two mechanisms such as the gas release because of the pressure drop (momentum transfer) and the condensation process (mass transfer) and their combination in different sequence and time. This “plateau” in the pressure time history values is formed at the remote locations due to increase of gas-to-liquid friction. The friction is increased because the liquid is propagated up towards the pressurized pipe section due to the condensation.
20
Page 20 of 38
Acknowledgement
The author would like to acknowledge the support of PETROSOFT-DC for giving
ip t
possibility to use the computer program (GDP code). The mathematical modeling and numerical analysis are performed using the GDP model. The author would like to acknowledge that
cr
permission was given by TransCanada (TCPL) to use the photos taken at the Didsbury Gas
an
us
Dynamic Test facility in Didsbury, Alberta (Canada).
M
References
Bendiksen, K.H., Maines, D., Moe, R., Nuland, S., 1991. The Dynamic Two-Fluid Model
te
d
OLGA: Theory and Application. SPE Production Engineering 6 (2), 171-180.
Blasius, P.R.H., 1913. Das Aehnlichkeitsgesetz bei Reibungsvorgangen in Fluessigkeiten.
Ac ce p
Forschungsheft. 131, 1-41.
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2003. Measurement of decompression wave speed in rich gas mixtures using a decompression tube. American Gas Association Proceedings -AGA-2003.
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2004. Determination of decompression wave speed in rich gas mixtures. The Canadian J. of Chemical Eng. 82, 880-891.
21
Page 21 of 38
Botros, K.K., Geerligs, J., Zhou, J., Glover, A., 2007. Measurements of flow parameters and decompression wave speed follow rapture of rich gas pipelines, and comparison with GASDECOM. Int. J. of Pressure Vessels and Piping. 84, 358-367.
ip t
Botros, K.K., Geerligs, J., Eiber, R.J., 2010a. Measurement of decompression wave speed in rich gas mixtures at high pressures (370 bars) using a specialized rupture tube. J. of Pressure Vessel
cr
Tech. 132, 051303-15.
us
Botros, K.K., Geerligs, J., Rothwell, B., Carlson, L., Fletcher, L., Venton, P., 2010b.
an
Transferability of decompression wave speed measured by a small-diameter shock tube to full size pipelines and implications for determining required fracture propagation resistance. Int. J. of
M
Pressure Vessels and Piping. 87, 681-695.
d
Burlutskiy, E., 2012a. Numerical analysis on rapid decompression in conventional dry gases
Ac ce p
Madrid, Spain. 63, 287-291.
te
using one-dimensional modeling. Proc. of Int. Conf. on Fluid Mechanics and Thermal Eng.
Burlutskiy, E., 2012b. Mathematical modeling of fracture propagation control. Proc. of ASME Pressure Vessels and Piping Conf., Toronto, Canada, PVP2012-78025, 1-8.
Burlutskiy, E., 2013a. Mathematical modelling on rapid decompression in base natural gas mixtures under rupturing. Chemical Eng. Research and Design, 91(1), 63-69.
Burlutskiy, E., 2013b. A one-dimensional mathematical model of multi-component fluid flow in pipes and its application to rapid decompression in dry natural gas mixtures. International Journal of Pressure Vessels and Piping, DOI: 10.1016/j.ijpvp.2013.01.003, to be published.
22
Page 22 of 38
Chen, J.R., Richardson, S.M., Saville, G., 1995a. Modelling of two-phase blowdown from pipelines - I. A hyperbolic model based on variational principles. Chem. Eng. Science. 50(4), 695-713.
ip t
Chen, J.R., Richardson, S.M., Saville, G., 1995b. Modelling of two-phase blowdown from pipelines - II. A simplified numerical method for multi-component mixtures. Chem. Eng.
cr
Science. 50(13), 2173-2187.
us
Eiber, R.J., Bubenik, T.A., Maxey, W.A., 1993. GASDECOM, computer code for the
an
calculation of gas decompression speed. Fracture control for natural gas pipelines, PRCI Report
M
No. L51691.
Eiber, R.J., Carlson, L., Leis, B., 2004. Fracture control requirements for gas transmission
te
d
pipelines, Proc. of the Fourth Int. Conf. on Pipeline Tech. 437.
Lee, A.L., Gonzalez, M.N., Eakin, B.E., 1966. The viscosity of natural gases, J. of Petroleum
Ac ce p
Tech. 997-1000.
Mahgerefteh, H., Oke, A., Atti, O., 2006a. Modelling outflow following rupture in pipeline networks, Chem. Eng. Science. 61(6), 1811-1818.
Mahgerefteh, H., Oke, A., Rykov, Y., 2006b. Efficient numerical solution for highly transient flows, Chem. Eng. Science. 61(15), 5049-5056.
Meulbroek, P., 2002. Equations of state in exploration. Organic Geochemistry. 33, 613-634.
Patankar, S., 1980. Numerical heat transfer and fluid flow, Hemisphere Publishing, New York. 23
Page 23 of 38
Soave, G., 1979. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Science. 27, 1197-1203.
ip t
Taitel, Y., Dukler, A.E., 1976. A model for predicting flow regime transitions in horizontal and
cr
near horizontal gas-liquid flow. AIChE J. 22, 47-55.
Tam, V.H.Y., Cowley, L.T., 1988. Consequences of pressurised LPG releases: The Isle of Grain
an
us
full scale experiments. Proc. of the GASTECH-88 Conf., Kuala-Lumpur, Malaysia.
M
Wallis, G. B., 1969. One-dimensional two-phase flows, McGraw Hill, New York.
CP G D pipe
hG hF
Ac ce p
Nomenclature
te
d
Zucrow, M.J., Hoffman, J.D., 1976. Gas dynamics, Wiley, New York.
heat capacity of the gas phase (J/kg/K) diameter of the pipe (m) enthalpy of the gas mixture (J/kg) enthalpy of the liquid fraction of the gas mixture (J/kg)
M
molar mass of the gas mixture (g/mol)
MC
condensation rate (kg/ m 2 /s)
N
number of components in the gas mixture
P
total pressure (MPa)
PC i
critical pressure value of the component i of gas mixture (MPa) 24
Page 24 of 38
critical condensation pressure when the mass transfer in the gas starts (MPa)
Pin
initial pressure in the shock tube before the decompression (MPa)
PrG
dimensionless complex (-)
q G F
term representing a heat exchange between the gas and the liquid film phase
R
universal gas constant ( J / mol / K )
R G Wall
gas mixture to wall friction term
R F Wall
liquid phase to wall friction term
R G F
gas core and liquid film friction term
S
cross-sectional area of the pipe ( m 2 )
T
temperature of the gas mixture (K)
TF
temperature of the liquid film phase (K)
t
time (s)
M
an
us
cr
ip t
PCOND
critical temperature value of the component i of gas mixture (K)
TCOND
critical condensation temperature when the mass transfer in the gas starts (K)
Tin
initial temperature in the shock tube before the decompression (K)
UG
velocity of the gas mixture (m/s)
te
Ac ce p
UF
d
TC i
velocity of the liquid film phase (m/s)
UG F
velocity of the gas to the liquid film mass transfer (m/s)
V
volume of the gas mixture ( m 3 )
Z
compressibility factor of the gas mixture
z
axial co-ordinate of the decompression tube (m)
zi
mole fraction of the component i of gas mixture
G
volume fraction of the gas mixture
F
volume fraction of the liquid film phase 25
Page 25 of 38
mass transfer term due to condensation from gas to liquid
internal surface roughness (m)
G
coefficient of gas thermal conductivity (W/m/K)
G
gas mixture viscosity ( Pa s )
F
liquid phase viscosity ( Pa s )
G Wall
gas to wall friction coefficient
F Wall
liquid film to wall friction coefficient
G F
gas to liquid film friction coefficient
perimeter of the pipe (m)
i
perimeter of the gas core at core-annular flow regime (m)
G
density of the gas mixture ( kg
G in
initial density of the gas mixture before the decompression ( kg
F
density of the liquid film phase ( kg
G Wall
gas mixture to wall friction term
F Wall
liquid film to wall friction term
G F
gas mixture to liquid film friction term
i
acentric factor of the component i of gas mixture
an
us
cr
ip t
G F
M
/m3)
/m3)
Ac ce p
te
d
/m3)
26
Page 26 of 38
Fig. 1 – Solving algorithm scheme.
Fig. 2 – Photo of the shock tube facility at TransCanada Gas Dynamic Test facility, Didsbury,
cr
Fig. 3 - Schematic of the decompression tube.
ip t
Alberta, Canada (Botros et. al., 2007).
an
selected PT locations.
us
Fig. 4 - Measured (Botros et. al., 2007) and predicted temperature time history values at the
M
Fig. 5 - Measured (Botros et. al., 2007) and predicted pressure time history values at the selected
d
PT locations.
te
Fig. 6 - Measured (Botros et. al., 2007) and predicted decompression wave speed as a function
Ac ce p
of pressure ratio.
Fig. 7 – Predicted propagation of the liquid phase front and gas decompression wave for the specified pressure values.
Fig. 8 - Predicted time history values of the slip between velocities of the liquid phase and the gas mixture for the selected tube elevations.
Fig. 9 - Predicted time history values of the liquid and gas velocities for the selected tube elevations.
Fig. 10 - Predicted pressure time history values for the selected tube elevations. 27
Page 27 of 38
Rupture disc
PT-2 / T-2
PT-3 / T-3
PT-4
PT-6
PT-8 / T-8
PT-10
0
0.06
0.24
0.44
0.84
1.64
4.04
cr
ip t
Distance (m)
Tin
N2
CO2
9.949
268.01
0.408
0.553
C1
C2
te
Pin
d
M
an
probe locations.
us
Table 1 – Distances between the rupture disc and Pressure Transducer (PT)/ Temperature (T)
21.406
i-C4
n-C4
i-C5
n-C5
C6+
9.080
0.026
0.014
0.002
0.002
0.001
Ac ce p
68.509
C3
Table 2 - Initial pressure (MPa), temperature (K) and gas composition (mole %).
28
Page 28 of 38
P Pin and T Tin
Solve:
SRK-EOS:
an
G f (Pin , Tin )
Film mass continuity Film mom. continuity Gas mass continuity Gas mom. continuity
M
Solve:
us
Input:
cr
ip t
Figure_1
-
Eq. (2) Eq. (4) Eq. (1) Eq. (3)
ed
Calculate: G f ( momentum , mass transfer ) SRK-EOS:
pt
Solve:
Ac
ce
Solve:
new old Pnew f (G ,T )
Film Enthalpy - Eq. (6) Gas Enthalpy - Eq. (5)
ti t i 1 t yes
t i t END no END
Page 29 of 38
an
us
cr
ip t
Figure_2
Ac
ce pt
ed
M
Fig. 2
Page 30 of 38
M
an
us
cr
ip t
Figure_3
Ac
ce pt
ed
Fig. 3
Page 31 of 38
Figure_4
0 PT-2 Measured PT-2 Predicted PT-8 Measured PT-8 Predicted
-5
-15 -20
ip t
-25 -30 -35 10
20
30 Time (ms)
40
50
an
us
-40
cr
o
Temperature ( C )
-10
Ac
ce pt
ed
M
Fig. 4
Page 32 of 38
Figure_5
10
PT-8 Measured PT-8 Predicted PT-10 Measured PT-10 Predicted
8
ip t
7 6 5 40 60 Time (ms)
80
an
us
20
cr
Pressure ( MPa )
9
Ac
ce pt
ed
M
Fig. 5
Page 33 of 38
Figure_6
1.0
0.8 0.7 0.6
Measured 2 Pred. (MC= 12.25 kg/m /s)
0.5
Pred. (MC= 49 kg/m /s)
ip t
Pressure Ratio ( - )
0.9
2
2
0
50 100 150 200 250 300 Decompression Wave Speed (m/s)
350
an
us
0.4
cr
Pred. (MC= 6.125 kg/m /s)
Ac
ce pt
ed
M
Fig. 6
Page 34 of 38
Figure_7
3.5
Liquid film Gas (9.8 MPa) Gas (9.0 MPa) Gas (8.0 MPa) Gas (7.0 MPa) Gas (6.0 MPa)
2.5 2.0
ip t
1.5 1.0 0.5 0
2
4
6 Time (ms)
8
10
an
us
0.0
cr
Pipe Axial Elevation ( m )
3.0
Ac
ce pt
ed
M
Fig. 7
Page 35 of 38
Figure_8
80 Predicted - z = 0.1 m Predicted - z = 0.5 m Predicted - z = 1.0 m
70
50 40 30
ip t
Slip Velocity ( m/s )
60
20 10
2
4 6 Time (ms)
8
10
an
us
0
cr
0
Ac
ce pt
ed
M
Fig. 8
Page 36 of 38
Figure_9
120
Gas ( z = 0.1 m) Liquid ( z = 0.1 m) Gas ( z = 0.5 m) Liquid ( z = 0.5 m) Gas ( z = 1.0 m) Liquid ( z = 1.0 m)
100
60 40
ip t
Velocity (m/s)
80
0 1
2
3
4 5 6 Time (ms)
7
8
9
10
an
us
0
cr
20
Ac
ce pt
ed
M
Fig. 9
Page 37 of 38
Figure_10
11 z = 0.1 m z = 0.5 m z = 1.0 m
10
8 7 6
ip t
Pressure (MPa)
9
5
0
1
2
3
4 5 6 Time (ms)
7
8
9
10
an
us
3
cr
4
Ac
ce pt
ed
M
Fig. 10
Page 38 of 38
S0957-5820(13)00014-1 http://dx.doi.org/doi:10.1016/j.psep.2013.03.001 PSEP 351
To appear in:
Process Safety and Environment Protection
Received date: Revised date: Accepted date:
21-11-2012 5-3-2013 12-3-2013
Please cite this article as: Burlutskiy, E., Numerical analysis of phase behavior during rapid decompression of rich natural gases, Process Safety and Environment Protection (2013), http://dx.doi.org/10.1016/j.psep.2013.03.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical analysis of phase behavior during rapid decompression of rich natural gases
Evgeniy Burlutskiy
ip t
A*STAR Institute of High Performance Computing (IHPC), 1 Fusionopolis Way, 138632 Singapore
cr
Tel: (65) 6419 1386 Fax: (65) 6467 4350 e-mail: [email protected]
us
Highlights 1D transient model of gas-liquid two-phase flow in a shock tube is developed
an
Approach to model the condensation in a compressible gas mixture is proposed The model is validated on the experimental data on RGD in rich natural gases
M
Physical reasons form the “condensation plateau” is investigated numerically
d
Abstract
te
The effect of the condensation process on the gas and liquid phase behavior during rapid decompression of rich natural gases is studied in the paper numerically. A one-dimensional
Ac ce p
mathematical model of transient thermal two-phase flow of compressible multi-component natural gas mixture and liquid phase in a shock tube is developed. The set of mass, momentum and enthalpy conservation equations is solved for the gas and liquid phases. The approach to model a liquid condensation process during rapid decompression of rich natural gas mixture is proposed. The mass transfer between the gas and the liquid is taken into account by introducing the appropriate terms into the governing equations. Thermo-physical properties of multicomponent natural gas mixture are calculated by solving the Equation of State (EOS) in the form of the Soave-Redlich-Kwong (SRK-EOS) model. The proposed liquid condensation model is integrated into the GDP model. A simple case of GDP model, where the liquid was not considered, was extensively validated on base and dry natural gases. The proposed two-phase
1
Page 1 of 38
model is validated against the experiments where the decompression wave speed was measured in rich natural gases at low temperature. It shows a good agreement with the experimental data.
KEY WORDS: Mathematical modeling; multiphase flow; Mass transfer; Condensation; Gases;
cr
ip t
Rapid decompression.
us
Introduction
an
Gas processing plants as well as basic fossil-fuel customers move further away from source production fields and carbon capture and storage (CCS) places. Rich and lean natural
M
gases and carbon dioxide (CO2) mixtures are transported by pipeline systems which are very long in length. Fractures in pipes caused by ruptures bring much more severe consequences
d
nowadays because of high operating pressure. A crack may easy propagate and destroy the
te
entire pipeline. The fracture propagation control in pipeline systems is usually based on the employment of the Battelle two-curve method (BTCM) developed by the Battelle Columbus
Ac ce p
Laboratories in order to determine the minimum of fracture arrest toughness (Eiber et. al., 1993, 2004). The crack propagation is arrested if the decompression wave speed in the gas is quicker than the fracture propagation velocity in the pipeline wall. Analytical models such as GASDECOM (Eiber et. al., 1993) are usually used by gas pipeline engineers in order to calculate the thermo-physical properties of natural gases. The program predicts the decompression wave speed in natural gases with a reasonably high level of accuracy (Botros et. al., 2003, 2007) using different equation of state (EOS) models. The effect of gas-to-wall friction is not taken into account because the fluid transport equations are not considered. The GASDECOM's predictive performance is poor if a relatively narrow shock tube (where the inner diameter is less than 5-10 cm) having a significant roughness of the internal surface is considered. The GASDECOM model does not work in the case of gas-liquid two-phase flows 2
Page 2 of 38
and if liquid condensation which is usually observed in the case of rich natural gases and carbon dioxide (CO2) mixtures is occurred. The effect of liquid on the decompression process may be accounted for in GASDECOM manually by introducing an average phase properties but the accuracy of such solutions is poor. However, more detailed information about time history
ip t
values and spatial distribution of basic flow parameters is also important. A mathematical model of fast transient multi-component fluid flows in a shock tube is highly desirable in this case in
cr
order to simulate a complex fluid dynamics. The problem became multiphase if the
us
condensation is occurred in the case of rich natural gases and carbon dioxide (CO2) mixtures.
an
Experimental investigation of single- and two-phase fluid dynamics in a shock tube is much easy nowadays because the accuracy of experimental instruments and the computer power
M
are significantly increased. First full scale experiment on the decompression (Tam and Cowley, 1988) were performed by British Petroleum and Shell Oil on the Isle of Grain (UK). A
d
horizontal pipe having a length of 100 meters and a diameter of 15 cm was used. The initial
te
pressure was varied from 0.8 MPa to 2.1 MPa in this study. The LPG mixture consisting of 95% mole propane and 5% mole butane was used as the fluid phase. Wide experimental
Ac ce p
measurements of the decompression wave speed in rich, base and dry natural gas mixtures were conducted by Trans Canada Pipe Lines (Botros et. al., 2003, 2004, 2007, 2010a, 2010b). Most of the measurements were made by using small-diameter shock tubes, where the friction influences on the flow behavior much stronger as compared to large-diameter pipes. The effect of shock tube diameter, gas mixture composition and fluid thermo-physical properties was examined experimentally. Pressure values were varied in those studies in the range between 10 MPa and 37 MPa (Botros et. al., 2003, 2010a). The temperature was varied from normal values to low (Botros et. al., 2003, 2007, 2010a). Time history data collections of basic flow parameters were obtained experimentally for rich natural gases in Botros et. al. (2007).
3
Page 3 of 38
Mathematical modeling of multi-phase fluid dynamics in a shock tube having natural gases at high pressure was also conducted. A simple method of solving the problem of multiphase blowdown flow from short pipes was proposed by Chen et. al. (1995a, 1995b). The set of two-phase flow governing equations was solved for multi-component fluid mixtures. The
ip t
proposed model was validated on the experiments of rapid decompression of simple gases. Numerical predictions showed a reasonable fit to the experiments of Tam and Cowley (1988). A
cr
numerical model of fluid dynamics following full-bore rupture in pipe systems having multicomponent fluid mixtures was proposed in Mahgerefteh et. al. (2006a, 2006b). The set of
us
governing equations was solved by using the method of characteristics (Zucrow and Hoffman,
an
1976). The model assumes a homogeneous equilibrium of the fluid flow, where the constituent phases are at thermal and mechanical equilibrium. The model focuses on simulation of complex
M
large-scale pipeline systems having bends, branches and couplings. The model was validated through the experiments of Tam and Cowley (1988) and a hypothetical pipeline system having a
d
length of 25 km successfully. A condensation process was not considered in those studies. A
te
one-dimensional code OLGA developed by SPT-group (Benediksen et. al., 1991) was examined in (Botros et. al., 2007) as a numerical solver for the problem of rapid decompression of dry and
Ac ce p
base natural gas mixtures. The attempt was un-successful and simulations showed a significant over-prediction of the pressure time history and a poor comparison (Burlutskiy 2013b) with the experimental data.
More general approach was applied in Burlutskiy (2012a, 2013a) for single gas and gas-
liquid flows in a shock tube (GDP model). The set of governing equations was solved for the gas and liquid phases together with taking into account of the Equation of State (EOS) model. The pressure time evolution in a shock tube for different axial elevations was predicted directly using the mass balance approach (Burlutskiy, 2013a, 2013b). The fracture propagation control was also successfully performed using the proposed GDP model (Burlutskiy, 2012b) in the frame of the Battelle two-curve method (BTCM), where two curves such as the crack velocity and the 4
Page 4 of 38
decompression wave speed were numerically simulated and compared with each other. Gasliquid system where the liquid is represented by the film phase and the condensation process were not considered in those works but the GDP model showed a great potential in modelling of
ip t
the liquid condensation process.
Numerical investigation of the condensation effect on the gas and liquid phase’s
cr
behavior during rapid decompression of rich natural gases is presented in the paper. The analysis is performed using a one-dimensional mathematical model of transient thermal two-phase flow
us
of compressible multi-component gas mixture and liquid film phase in a shock tube. The
an
approach to model the process of liquid condensation during rapid decompression of rich natural gases as well as the numerical algorithm of solving the problem is proposed in the paper. The
M
proposed model is validated against the experiments (Botros et. al., 2007) where the
te
d
decompression wave speed in rich natural gases was measured.
Ac ce p
One-dimensional mathematical model of gas-liquid two-phase flow in a shock tube
A one-dimensional transient mathematical model of compressible thermal gas-liquid
two-phase flow in a shock tube is developed. The set of mass, momentum and enthalpy conservation equations is solved for the gas and liquid phases. The liquid phase is assumed to be in the form of liquid film because the internal pipe diameter in the considered case is very small (49.325 mm) as compared to industrial pipeline diameters (about 1 m). The mass transfer has a form of fast condensation process from the gaseous phase into the liquid film phase. The natural gas mixture is considered as the gas core phase in this system. The set of governing equations in general form is written as (Wallis, 1969):
5
Page 5 of 38
G G G G U G G F t z
(1)
F F F F U F t z
(2)
G F
2 G G U G G G U G P G t z z
F
P z
R GF
G F U G F
R G F
R G Wall
R F Wall
(3)
cr
F F U F F F U 2F t z
G F U G F
ip t
an
P (U UG )2 P q G F G F G F F UF z 2 t R G F ( U F U G ) R F Wall U F
(5)
(6)
M
F F h F F F U Fh F t z
us
P G G h G G G U G h G P (U U G )2 q G F G F G F G UG t z z 2 t R G F ( U F U G ) R G Wall U G
(4)
Here, G , F are the volume fraction of gas mixture and liquid film, respectively; G ,
d
F are the density of the gas and the film; UG , UF are the velocity of the gas and the film,
te
respectively; P is the total pressure; R G Wall , R F Wall are friction terms, which represent the
Ac ce p
gas-to-wall and the film-to-wall friction, respectively; R G F is the term, which represents the friction between the gas and film phases; h G , h F are the enthalpy of the gas and the film; q G F is the term, which represents a heat exchange between the gas phase and the liquid film phase; G F is the term, which represents a condensation from the gaseous phase into the liquid phase;
t is the time; z is the axial co-ordinate. Last terms in the enthalpy conservation equations (5)-(6) represent additional contribution into the system’s enthalpy due to the friction force. It is assumed that the velocity of gas fraction, which leaves the gaseous phase during the condensation process and passes into the liquid state, is equal to the velocity of the liquid fraction (i.e. U G F U F ) during the condensation.
6
Page 6 of 38
The mass transfer term which describes a liquid condensation process is represented by the following correlation:
i MC S
(7)
ip t
G F
Here, M C is the condensation rate. The mass transfer term (7) is a universal correlation which is
cr
dependent from the shock tube inner diameter (i.e. i D pipe ; S D 2pipe / 4 ); i is the
us
perimeter of the gas core phase at core-annular flow regime (i.e. the gas core-to-liquid film interface); S is the cross-sectional area of the pipe. Therefore, it can be easy extrapolated into
an
other cases. The condensation rate ( M C ) is a constant in the proposed study and its value will be determined numerically through a better fit to the experimental data. A dependence of the
M
condensation rate ( M C ) from the various flow parameters such as the gas mixture composition, temperature, pressure and others will be investigated in other studies. The friction term, which
G Wall S
(8)
Ac ce p
R G Wall
te
d
represents gas-to-wall friction, is written in the Blasius form (Blasius, 1913):
G Wall
2 , Re G 1600 G Wall 64 / Re G G Wall G U G , 0.25 8 G Wall 0.316 / Re G , Re G 1600
ReG G UG DpipeG / G
(9)
(10)
Here, is the perimeter of the pipe; G Wall is the gas-to-wall friction term; G Wall is
the friction coefficient, D pipe is the diameter of the pipe; G is the gas viscosity. The term, which describes the friction between the film and the wall, is written in the Taitel and Dukler form (Taitel and Dukler, 1976):
7
Page 7 of 38
R F Wall
F Wall S
F Wall F Wall F U 2F
(11) , Re F 1600 F Wall 64 / Re F F Wall 0.046 / Re 0F.2 , Re F 1600
,
(12)
ReF F UFDpipeF / F
ip t
(13)
cr
Here, F Wall is the liquid film-to-wall friction term; F Wall is the friction coefficient,
F is the viscosity of the liquid. The friction between the gas and the film is also written in the
us
Blasius form (Blasius, 1913):
i G F S
G F
, Re G 1600 G F 64 / Re G G F G U G U F ( U G U F ) , 0 . 25 8 G F 0.316 / Re G , Re G 1600
(14)
(15)
d
M
an
R GF
te
Here, G F is the gas-to-film friction term; G F is the friction coefficient. The heat
Ac ce p
transfer between the gas phase and the liquid film phase is written as (Wallis, 1969):
qGF
i i G Nu G F (T TF ) S 4 S(1 F )
(16)
0.5 0.33 is the Here, CP G is the heat capacity of the gas phase; NuG F 2 0.6 ReG F PrG
Nusselt number of the flow; G is the coefficient of heat conductivity of the gas phase.
Thermo-physical properties of natural gas mixtures
8
Page 8 of 38
Thermo-physical fluid properties are modeled by solving of the Equations of State (EOS) in the form of the Soave-Redlich-Kwong (SRK-EOS) model (Soave, 1979). It has been chosen due to better performance as compared to all other EOS models. The set of equations and
a V V b
N N
a zi z j a i a j 1 kij i 1j1
(17)
;
N
b z i bi i 1
2
an
R 2 TC2 i T 1 m i 1 a i 0.42748 PC i TC i
RTC i PC i
M
bi 0.08664
cr
V b
us
RT
P
ip t
correlations (SRK-EOS) may be written as (Soave, 1979):
(19)
(20)
(21)
d
m i 0 . 48 1 . 574 i 0 . 176 i2
(18)
te
Here, V is the volume of the gas mixture; N is the number of components in the gas
Ac ce p
mixture; T is the temperature of the gas mixture; R is the universal gas constant; i is the acentric factor of the component i; PC i , TC i are critical values of the pressure and temperature, correspondently; zi is the mole fraction of the component i. The compressibility factor (Z) of gas mixture is calculated from the following equation (Soave, 1979):
Z3 Z2 A B B2 Z AB 0 A
aP R 2T 2
P G
;
B
(22)
bP RT
(23)
R TZ M
(24)
9
Page 9 of 38
Here, M is the molar mass of a substance. The viscosity of gas mixture is calculated by using the Lee-Gonzales-Eakin (LGE) correlation (Lee et. al., 1966).
ip t
Solving algorithm of two-phase problem
cr
The algorithm of solving of the set of one-dimensional transient governing equations of the fluid mixture flow in a pipe is based on the Tri-Diagonal Matrix Algorithm (TDMA), also
us
known as the Thomas algorithm (Patankar, 1980). It is a simplified form of Gaussian
an
elimination that can be used to solve tri-diagonal systems of equations. The computational pipe is divided into N parts, which are called the control volumes. The set of unsteady governing
M
equations is transformed into the standard form of the discrete analog of the tri-diagonal system (Patankar, 1980) by using the fully implicit numerical scheme and it is calculated at each control
d
volume of the computational pipe. The equation is reduces to the steady state discretization
te
equation if the time step goes to infinity.
Ac ce p
The gas and liquid initial velocities are set to zero ( U G i = 0 m/s; U L i = 0 m/s; i=1, ..., N;
N is the number of control volumes of the pipe) in each control volume of the computational pipe before the decompression. The pressure and the temperature of the gas mixture are set to specified initial values. The pressure outside the tube is set to the atmospheric ( PN Pin ; PN 1 = 0.1 MPa; PN >> PN 1 ; N+1 is a ghost cell). The initial volume fraction values of the liquid are set to zero. The initial density of the gas mixture is calculated by using the SRK-EOS model. This is the direct way of using the SRK-EOS model, where the initial density of the gas mixture is a function of the initial pressure and temperature (fig. 1). The following assumptions are made in order to model the compressible nature of the process:
10
Page 10 of 38
The gas mixture density is assumed to be constant when the mass and momentum conservation equations are solved for both phases. A change in the gas density is determined after when the set of governing equations is solved for the considered time step.
The time step of
2 10 5
ip t
The liquid phase density is a constant during the decompression process in the study. sec is selected for simulations in order to provide a continuous
us
cr
change of the gas density.
Main iteration loop starts from the momentum transfer predictions (fig. 1). Velocity
an
values of the gas and liquid film phases are calculated by solving the momentum conservation equations of each phase. The mass continuity equations are solved for both phases after solving
M
the momentum transfer. The momentum and mass continuity equations are solved for the liquid film phase if the liquid is already occurred. The volume fraction of the gas is conceived as a
d
local variable for one time step. It is set to 1 (one) at each control volume of the pipe before
te
solving the momentum and mass continuity equations. When the time step is calculated, a new density of the gas mixture is determined. The old gas density values from the previous time step
Ac ce p
are multiplied on the gas volume fraction values from new time step in order to account for the momentum change for the considered time step. Further, the gas volume fraction values are set to 1 (one) at each control volume for a new time step calculations again. The condensation process is taken into account by introducing the appropriate terms into the governing equations. The mass balance between the gas and the liquid phase for the selected control volume is considered, keeping the total mass of both phases in the selected control volume constant for the considered time step.
When all changes in the gas density is taken into account, new pressure values are determined thereafter. These values are calculated by using the SRK-EOS model determining the inverse correspondence, which may be written schematically as (fig. 1): 11
Page 11 of 38
new old Pnew f (G ,T )
(25)
So, new pressure values are determined at each control volume of the pipe. New
ip t
temperature values for both phases are calculated thereafter by solving the enthalpy conservation equations for the gas and liquid phases. Thus, the density, pressure, temperature and velocity
cr
values are new. This procedure is repeated again for a new time step if the total calculation time
us
is smaller than selected prediction time. The mass transfer is started if the pressure and temperature values of the gas mixture are lower than the critical condensation values (i.e.
an
Pi PCOND and Ti TCOND ) which are dependent from gas mixture composition. The condition
of liquid occurring in the shock tube is checked at every time step before solving the set of
M
governing equations. The proposed model was developed under the research project “Fast Transient Flows of Multi-Component Natural Gas Mixture and Liquid in Pipelines” in
d
PETROSOFT-DC. It was implemented into the FORTRAN computer code named the “Gas
Ac ce p
te
Decompression Program” (GDP code).
Numerical simulation of rapid decompression of rich natural gas mixture
A simple case of the proposed GDP model, where the condensation was not considered,
has been validated against the experimental data on rapid decompression of a shock tube having dry and base natural gases (Botros et. al., 2003) and presented in Burlutskiy (2013a, 2013b). Single-phase flows of natural gases at thermal equilibrium were considered. The proposed twophase model, where the condensation process is taken into account and integrated into the GDP code, is validated against the decompression wave speed measurements were made by Trans Canada Pipe Lines (TCPL) at TCPL Gas Dynamic Test Facility in Didsbury, Alberta, Canada (Botros et. al., 2007) using the same experimental facility as in the case of dry and base natural 12
Page 12 of 38
gases (Burlutskiy, 2013a, 2013b). Fig. 2 shows a photo of the real experimental shock tube facility. The main test section is a shock tube, which has a total length of 172 m. The internal pipe diameter is 49.325 mm. It has a roughness better than 1 micro-meter (i.e. less than 1.0 10 6 6
m). Considering the roughness-to-diameter ratio ( / Dpipe 20.3 10 ;
1.0 10 6
m;
ip t
Dpipe 49.325103 m), the internal surface of the shock tube is determined as the “smooth” tube
cr
surface (Botros et. al., 2010b). The classical Blasius correlation (8)-(10) for the smooth pipes is allowed in this case. A rupture disc was placed at one end of the experimental shock tube, which
us
was upon rupturing. A decompression wave propagates up into the pressurized test section. High frequency responses Pressure Transducers (PT) were mounted into the tube part, which was
an
most close to the rupture end of the shock tube, in order to capture the time history of the expansion fan (Botros et. al., 2007). Distances between the pressure transducers and the rupture
M
disc are shown in table 1. The schematic of the decompression tube is shown in Fig. 3. Four
te
wave speed in those experiments.
d
pressure-time tracers (PT-2, PT-3, PT-4 and PT-6) were used to determinate the decompression
Ac ce p
Fig. 2 and Fig. 3 Table 1, 2
The computational decompression pipe has a length of 15 m. The set of governing
equations is solved on the mesh having 1000 grid nodes on the axial direction. The inner pipe diameter is 49.325 mm. One end of the shock tube is selected to be the closed end. The rupture disc condition is modeled in the other end of the pipe. The initial pressure ( Pin ), temperature ( Tin ), and natural gas mixture composition are shown in table 2 (Botros et. al., 2007). The considered rich natural gas mixture is in the gaseous state if the pressure is in the range between 8 and 10 MPa, and the temperature is between -5 o C and -10 o C . The gas mixture is expanded towards the open end of the shock tube if the rupture disc is opened. The pressure and the 13
Page 13 of 38
temperature of the gas are dropped down due to gas expansion. The mixture which is in the gaseous state at the beginning goes to two-phase region of the pressure-temperature envelope if the pressure is less than 8 MPa and the temperature is less than -10 o C . Below those critical P and T values, the mixture is transformed from single phase system (gas) to two-phase system
ip t
(gas-liquid). Those critical pressure and temperature parameters are dependent from the natural gas mixture composition. Thermo-physical properties of the considered two-phase gas-liquid
cr
system are calculated by employing the SRK-EOS model together with the flash algorithm
us
(Meulbroek 2002) in the conditions of two-phase region of pressure-temperature envelope. The liquid phase density is determined numerically using the SRK-EOS model and the flash kg / m 3 .
an
algorithm and it is equal to 350
M
Fig. 4
d
A collection of measured (Botros et. al., 2007) and predicted temperature time history
te
values at PT-2 and PT-8 locations are shown in fig. 4. The temperature was not collected at other considered PT’s (PT-8 and PT-10) at the experimental measurements. Experimental points
Ac ce p
are shown as symbols here. Broken curves represent calculated values. Predicted temperature time history values are in a good agreement with the measurements in most part of the pressure range. A decrease rate of predicted time-temperature curves is stronger as compared to measured points for pipe elevations near the rupture end (PT-2). The slope of the predicted curve at PT-8 pipe location is much better correlate with the experimental data. Hence, predicted decompressed temperature values (i.e. the level after the decompression) are in a good correlation with measured temperature minimum for every of the considered PT locations. The proposed model simulates the temperature field of the gas mixture with high order of agreement with measured values.
Fig. 5 14
Page 14 of 38
A comparison between measured (Botros et. al., 2007) and predicted pressure time history values at PT-8 and PT-10 locations is shown in fig. 5. Experimental points are shown as symbols. A continues line and broken curves represent calculated values. Measured pressure
ip t
values are decreased at each PT location of the shock tube if the decompression in the gas mixture is started at those PT locations. A decrease of the pressure is more rapid if the
cr
considered PT location closer to the rupture end of the pipe, and is slower for the PT’s far away from the rupture end due to strong effect of gas-to-wall friction. Hence, predicted pressure time
us
history values are in a good agreement with the measurements in most part of the pressure range.
an
The slope of the calculated time-pressure curves is in a good correlation with the experimental values. Moreover, predicted pressure values at the considered PT locations reach measured
M
pressure minimum after the decompression successfully.
te
d
Fig. 6
A correlation between the decompression wave speed and the pressure has to be
Ac ce p
determined in order to perform the Battelle two-curve analysis. The decompression wave speed is calculated by using predicted pressure time history values at PT-2, PT-3, PT-4 and PT-6 locations. Fig. 6 shows measured (Botros et. al., 2007) and predicted decompression wave speed values as a function of the pressure. Pressure values are normalized on the initial pressure in the shock tube before the decompression ( P / Pin ). The initial pressure is equal to 9.949 MPa in this study. Symbols correspond to measured values (Botros et. al., 2007). Curves are numerical simulations performed by using the proposed model, where the liquid condensation is taken into account. The condensation rate is determined numerically through a better fit to the experiments where the decompression wave speed in rich natural gas mixture was measured. The mass transfer effects on the decompression wave speed when the condensation process is occurred (i.e. the pressure ratio between 0.7 and 0.8 in fig. 6). Different values of the condensation rate 15
Page 15 of 38
( M C ) are examined numerically. Fig. 6 shows predictions of the decompression wave speed if the mass transfer rate ( M C ) is equal to 6.125, 12.25 and 49 (kg/m2/s). The best fit to measured values is observed (Fig. 6) if the condensation rate value is equal to 12.25 (kg/m2/s). Real values of the condensation rate were not measured experimentally. Therefore, the condensation rate
ip t
( M C ) equal to 12.25 (kg/m2/s) is chosen to be used in this study because it provides the best fit to measurements. A rapid decrease in the decompression wave speed is not predicted by the
cr
model if the mass transfer between phases is not taken into account or the condensation rate is
an
Fig. 7
us
too small.
M
Fig. 7 shows a development of the liquid phase front propagation up into the pressurized shock tube section predicted by the proposed model. It is compared with predicted lines which
d
represent the propagation of gas decompression wave front for few specified pressure values.
te
The pipe axial elevation in this figure shows the distance from the rupture end of the shock tube. The decompression wave in rich natural gas mixture propagates almost linearly if the pressure is
Ac ce p
above the critical condensation values (i.e. 9.8 MPa; 9.0 MPa; 8.0 MPa lines in fig. 7). The liquid film propagation shows a non-linear behavior in the beginning. The propagation velocity of the liquid phase front is significantly decreased during first 5-6 milliseconds as compared to the beginning. The liquid phase effects on the development of gas decompression wave propagation if pressure values below the critical (i.e. 7.0 MPa; 6.0 MPa lines in fig. 7) during first 2-3 milliseconds due to phases re-distribution and interaction between gas and liquid phases. Such non-linearity in the decompression wave propagation in rich natural gas mixture cannot be predicted directly by using of analytical models.
Fig. 8
16
Page 16 of 38
A flat plateau is observed in both pressure time history curves (fig. 5). This “condensation plateau” is observed if the condensation process in natural gases is occurred only and it is investigated in more details below. For that reason a few more figures which contain time history values at three different specified tube elevations are considered. Distances between
ip t
those elevations and the rupture disc are 0.1 m, 0.5 m and 1 m. Predicted time history values of a slip between velocities of the liquid film and the gas phase at those elevations are shown in fig. 8
cr
and analyzed. The slip is increased at every of those locations because the gas is accelerated when it is released. It consists from the gas only at the beginning because the liquid is occurred
us
later. The slip is increased until a certain level (15-20 m/s) and has a local maximum here. The
an
slip is decreased thereafter because the liquid is condensed and accelerated by the gas, and has a local minimum at the end of the condensation process (i.e. t = 2 ms for z = 0.1 m; t = 4 ms for = 0.5 m; t = 8 ms for z = 1 m). The condensation front keeps propagating up into the
M
z
pressurized pipe section. The decompression continues at those locations due to gas release.
d
However, the gas mixture loses much more kinetic energy than before by accelerating the liquid
te
phase. For the most remote locations ( z = 1 m) gas-to-liquid friction effects much stronger
Ac ce p
(fig. 8) because of a bigger distance from the rupture end, which is occupied by the liquid.
Fig. 9 and Fig. 10
Fig. 9 shows predicted velocity values of the gas and liquid phases which are collected at
the same pipe elevations (0.1 m, 0.5 m and 1 m). Fig. 10 shows pressure time evolution of the gas mixture at those elevations also. Both figures show that the “condensation plateau” in the pressure time history (fig 10, 7 MPa < P < 8 MPa) is bigger for the most distant axial elevation ( z = 1 m). The velocity of the gas mixture is accelerated from 0 m/s to 20 m/s for the elevation of 1 m due to gas expansion at the beginning (fig. 9,
z
= 1 m, 3 ms < t < 5 ms). A pressure
decrease for the considered pipe elevation of 1 m is slower as compared to the elevation of 0.1 m 17
Page 17 of 38
because the friction effects much stronger due to longer distance from the rupture end. A condensation is not occurred yet during this time period at the elevation of 1 m. The liquid condensation front is propagated up from the elevation of 0.6 m to the elevation of 0.8 m during this time period (fig. 7). The velocity of the gas does not show the acceleration and is almost
ip t
constant during the time between 4.5 and 7 ms. Gas expansion in this case is extremely slow due to strong increase of the distance where the liquid condensation is occurred already, and gas-to-
cr
liquid friction contribution, respectively. For the pipe elevation of 0.5 m the effect of the “condensation plateau” is weaker because the distance from the rupture end is shorter. For the
us
pipe elevation of 0.1 m the effect of the “condensation plateau” is absent because the friction is
an
very small in this case. Therefore, high friction between the gas and the liquid phase (gas-toliquid friction) forms the “condensation plateau” in the pressure. The condensation process is
M
started thereafter at the elevation of 1 m when the pressure and the temperature are dropped down and the critical pressure and temperature values are reached (fig. 9,
z
= 1 m, liquid, t = 7
z
= 1 m, liquid, 7 ms < t < 8.5 ms). A pressure decrease is weaker
te
mechanisms act both (fig. 9,
d
ms). The pressure is decreased again because the condensation and momentum transfer
again starting from the time of 8.5 ms because the only momentum transfer (gas release) acts on
Ac ce p
the gas decompression here. Fig. 9 and 10 show together that an increase in gas-to-liquid friction downstream of the considered pipe elevation, due to liquid front propagation toward the pressurized pipe section, forms the “condensation plateau” in the pressure time history.
The decompression in rich natural gas mixtures below the critical condensation
parameters is determined by the action of two mechanisms such as the gas release due to the pressure drop (the momentum transfer) and the liquid condensation process (the mass transfer). A rapid decrease in the decompression wave speed at condensation is dependent from the “condensation plateau” in the pressure time history, which is used to determined those values. Numerical analysis shows that the “condensation plateau” in the pressure time history values is formed due to the following reasons: 18
Page 18 of 38
First, the gas mixture leaves the considered pipe elevation and it is accelerated itself due to gas release. The liquid phase is occurred in the downstream due to the condensation and propagates up into the pressurized pipe section starting from the rupture end. The gas
ip t
interacts with the liquid in the downstream by losing the kinetic energy because it accelerates the liquid phase in the downstream. Such decrease of the gas kinetic energy
cr
in the downstream due to increase of gas-to-liquid friction (the momentum transfer) slows down of the gas acceleration by keeping the gas velocity constant. It reduces the
an
“plateau” in the pressure time history values.
us
pressure drop between the considered pipe elevation and the downstream, and forms the
M
At the later stage, the pressure and the temperature are dropped down enough and the critical pressure and temperature values are reached. The liquid phase is condensed in the
d
considered elevation and the decompression in this location is determined by the action
te
of two mechanisms such as the gas release (momentum transfer) and the condensation process (mass transfer) at the same time together. Therefore, a pressure decrease is
Ac ce p
significantly quicker at the considered elevation as compared to usual decrease rate caused by the momentum transfer only. The pressure decrease is slow again when the condensation is finished.
The effect of the liquid phase on the gas mixture flow behavior is much stronger if small-
diameter pipes and low temperatures are considered. The volume fraction of the liquid phase is much bigger in this case and the friction force effects on the flow behavior much stronger.
Conclusions
19
Page 19 of 38
The effect of the condensation process on the gas and liquid phase behavior during rapid decompression of rich natural gases is analyzed numerically. A one-dimensional mathematical model of transient thermal two-phase flows of compressible multi-component gas mixture and liquid phase in a shock tube is developed. The set of mass, momentum and enthalpy
ip t
conservation equations is solved for the gas and liquid phases. The approach to model the process of liquid condensation during rapid decompression of a compressible rich natural gas
cr
mixture as well as the numerical algorithm of solving the problem is proposed in the paper. The mass transfer process between the gas and the liquid is taken into account by introducing the
us
appropriate terms into the set of governing equations. The proposed model is validated on the
an
experiments where the decompression wave speed was measured in a shock tube having rich
M
natural gases at low temperature and it shows a good agreement with the experimental data.
Numerical analysis shows a non-linear development of the liquid phase front
d
propagation. The propagation speed of liquid phase front is significantly decreased during first
te
few milliseconds as compared to the propagation start. The liquid phase influences on the development of gas decompression wave propagation if the pressure values below the critical
Ac ce p
condensation values during first few milliseconds due to phases re-distribution and frictional interaction between the gas and liquid phases. Such non-linearity in the decompression wave development in natural gas mixtures cannot be predicted directly by using of analytical models.
Numerical analysis shows that the “condensation plateau” in the decompression wave
speed is taken place due to action of two mechanisms such as the gas release because of the pressure drop (momentum transfer) and the condensation process (mass transfer) and their combination in different sequence and time. This “plateau” in the pressure time history values is formed at the remote locations due to increase of gas-to-liquid friction. The friction is increased because the liquid is propagated up towards the pressurized pipe section due to the condensation.
20
Page 20 of 38
Acknowledgement
The author would like to acknowledge the support of PETROSOFT-DC for giving
ip t
possibility to use the computer program (GDP code). The mathematical modeling and numerical analysis are performed using the GDP model. The author would like to acknowledge that
cr
permission was given by TransCanada (TCPL) to use the photos taken at the Didsbury Gas
an
us
Dynamic Test facility in Didsbury, Alberta (Canada).
M
References
Bendiksen, K.H., Maines, D., Moe, R., Nuland, S., 1991. The Dynamic Two-Fluid Model
te
d
OLGA: Theory and Application. SPE Production Engineering 6 (2), 171-180.
Blasius, P.R.H., 1913. Das Aehnlichkeitsgesetz bei Reibungsvorgangen in Fluessigkeiten.
Ac ce p
Forschungsheft. 131, 1-41.
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2003. Measurement of decompression wave speed in rich gas mixtures using a decompression tube. American Gas Association Proceedings -AGA-2003.
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2004. Determination of decompression wave speed in rich gas mixtures. The Canadian J. of Chemical Eng. 82, 880-891.
21
Page 21 of 38
Botros, K.K., Geerligs, J., Zhou, J., Glover, A., 2007. Measurements of flow parameters and decompression wave speed follow rapture of rich gas pipelines, and comparison with GASDECOM. Int. J. of Pressure Vessels and Piping. 84, 358-367.
ip t
Botros, K.K., Geerligs, J., Eiber, R.J., 2010a. Measurement of decompression wave speed in rich gas mixtures at high pressures (370 bars) using a specialized rupture tube. J. of Pressure Vessel
cr
Tech. 132, 051303-15.
us
Botros, K.K., Geerligs, J., Rothwell, B., Carlson, L., Fletcher, L., Venton, P., 2010b.
an
Transferability of decompression wave speed measured by a small-diameter shock tube to full size pipelines and implications for determining required fracture propagation resistance. Int. J. of
M
Pressure Vessels and Piping. 87, 681-695.
d
Burlutskiy, E., 2012a. Numerical analysis on rapid decompression in conventional dry gases
Ac ce p
Madrid, Spain. 63, 287-291.
te
using one-dimensional modeling. Proc. of Int. Conf. on Fluid Mechanics and Thermal Eng.
Burlutskiy, E., 2012b. Mathematical modeling of fracture propagation control. Proc. of ASME Pressure Vessels and Piping Conf., Toronto, Canada, PVP2012-78025, 1-8.
Burlutskiy, E., 2013a. Mathematical modelling on rapid decompression in base natural gas mixtures under rupturing. Chemical Eng. Research and Design, 91(1), 63-69.
Burlutskiy, E., 2013b. A one-dimensional mathematical model of multi-component fluid flow in pipes and its application to rapid decompression in dry natural gas mixtures. International Journal of Pressure Vessels and Piping, DOI: 10.1016/j.ijpvp.2013.01.003, to be published.
22
Page 22 of 38
Chen, J.R., Richardson, S.M., Saville, G., 1995a. Modelling of two-phase blowdown from pipelines - I. A hyperbolic model based on variational principles. Chem. Eng. Science. 50(4), 695-713.
ip t
Chen, J.R., Richardson, S.M., Saville, G., 1995b. Modelling of two-phase blowdown from pipelines - II. A simplified numerical method for multi-component mixtures. Chem. Eng.
cr
Science. 50(13), 2173-2187.
us
Eiber, R.J., Bubenik, T.A., Maxey, W.A., 1993. GASDECOM, computer code for the
an
calculation of gas decompression speed. Fracture control for natural gas pipelines, PRCI Report
M
No. L51691.
Eiber, R.J., Carlson, L., Leis, B., 2004. Fracture control requirements for gas transmission
te
d
pipelines, Proc. of the Fourth Int. Conf. on Pipeline Tech. 437.
Lee, A.L., Gonzalez, M.N., Eakin, B.E., 1966. The viscosity of natural gases, J. of Petroleum
Ac ce p
Tech. 997-1000.
Mahgerefteh, H., Oke, A., Atti, O., 2006a. Modelling outflow following rupture in pipeline networks, Chem. Eng. Science. 61(6), 1811-1818.
Mahgerefteh, H., Oke, A., Rykov, Y., 2006b. Efficient numerical solution for highly transient flows, Chem. Eng. Science. 61(15), 5049-5056.
Meulbroek, P., 2002. Equations of state in exploration. Organic Geochemistry. 33, 613-634.
Patankar, S., 1980. Numerical heat transfer and fluid flow, Hemisphere Publishing, New York. 23
Page 23 of 38
Soave, G., 1979. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Science. 27, 1197-1203.
ip t
Taitel, Y., Dukler, A.E., 1976. A model for predicting flow regime transitions in horizontal and
cr
near horizontal gas-liquid flow. AIChE J. 22, 47-55.
Tam, V.H.Y., Cowley, L.T., 1988. Consequences of pressurised LPG releases: The Isle of Grain
an
us
full scale experiments. Proc. of the GASTECH-88 Conf., Kuala-Lumpur, Malaysia.
M
Wallis, G. B., 1969. One-dimensional two-phase flows, McGraw Hill, New York.
CP G D pipe
hG hF
Ac ce p
Nomenclature
te
d
Zucrow, M.J., Hoffman, J.D., 1976. Gas dynamics, Wiley, New York.
heat capacity of the gas phase (J/kg/K) diameter of the pipe (m) enthalpy of the gas mixture (J/kg) enthalpy of the liquid fraction of the gas mixture (J/kg)
M
molar mass of the gas mixture (g/mol)
MC
condensation rate (kg/ m 2 /s)
N
number of components in the gas mixture
P
total pressure (MPa)
PC i
critical pressure value of the component i of gas mixture (MPa) 24
Page 24 of 38
critical condensation pressure when the mass transfer in the gas starts (MPa)
Pin
initial pressure in the shock tube before the decompression (MPa)
PrG
dimensionless complex (-)
q G F
term representing a heat exchange between the gas and the liquid film phase
R
universal gas constant ( J / mol / K )
R G Wall
gas mixture to wall friction term
R F Wall
liquid phase to wall friction term
R G F
gas core and liquid film friction term
S
cross-sectional area of the pipe ( m 2 )
T
temperature of the gas mixture (K)
TF
temperature of the liquid film phase (K)
t
time (s)
M
an
us
cr
ip t
PCOND
critical temperature value of the component i of gas mixture (K)
TCOND
critical condensation temperature when the mass transfer in the gas starts (K)
Tin
initial temperature in the shock tube before the decompression (K)
UG
velocity of the gas mixture (m/s)
te
Ac ce p
UF
d
TC i
velocity of the liquid film phase (m/s)
UG F
velocity of the gas to the liquid film mass transfer (m/s)
V
volume of the gas mixture ( m 3 )
Z
compressibility factor of the gas mixture
z
axial co-ordinate of the decompression tube (m)
zi
mole fraction of the component i of gas mixture
G
volume fraction of the gas mixture
F
volume fraction of the liquid film phase 25
Page 25 of 38
mass transfer term due to condensation from gas to liquid
internal surface roughness (m)
G
coefficient of gas thermal conductivity (W/m/K)
G
gas mixture viscosity ( Pa s )
F
liquid phase viscosity ( Pa s )
G Wall
gas to wall friction coefficient
F Wall
liquid film to wall friction coefficient
G F
gas to liquid film friction coefficient
perimeter of the pipe (m)
i
perimeter of the gas core at core-annular flow regime (m)
G
density of the gas mixture ( kg
G in
initial density of the gas mixture before the decompression ( kg
F
density of the liquid film phase ( kg
G Wall
gas mixture to wall friction term
F Wall
liquid film to wall friction term
G F
gas mixture to liquid film friction term
i
acentric factor of the component i of gas mixture
an
us
cr
ip t
G F
M
/m3)
/m3)
Ac ce p
te
d
/m3)
26
Page 26 of 38
Fig. 1 – Solving algorithm scheme.
Fig. 2 – Photo of the shock tube facility at TransCanada Gas Dynamic Test facility, Didsbury,
cr
Fig. 3 - Schematic of the decompression tube.
ip t
Alberta, Canada (Botros et. al., 2007).
an
selected PT locations.
us
Fig. 4 - Measured (Botros et. al., 2007) and predicted temperature time history values at the
M
Fig. 5 - Measured (Botros et. al., 2007) and predicted pressure time history values at the selected
d
PT locations.
te
Fig. 6 - Measured (Botros et. al., 2007) and predicted decompression wave speed as a function
Ac ce p
of pressure ratio.
Fig. 7 – Predicted propagation of the liquid phase front and gas decompression wave for the specified pressure values.
Fig. 8 - Predicted time history values of the slip between velocities of the liquid phase and the gas mixture for the selected tube elevations.
Fig. 9 - Predicted time history values of the liquid and gas velocities for the selected tube elevations.
Fig. 10 - Predicted pressure time history values for the selected tube elevations. 27
Page 27 of 38
Rupture disc
PT-2 / T-2
PT-3 / T-3
PT-4
PT-6
PT-8 / T-8
PT-10
0
0.06
0.24
0.44
0.84
1.64
4.04
cr
ip t
Distance (m)
Tin
N2
CO2
9.949
268.01
0.408
0.553
C1
C2
te
Pin
d
M
an
probe locations.
us
Table 1 – Distances between the rupture disc and Pressure Transducer (PT)/ Temperature (T)
21.406
i-C4
n-C4
i-C5
n-C5
C6+
9.080
0.026
0.014
0.002
0.002
0.001
Ac ce p
68.509
C3
Table 2 - Initial pressure (MPa), temperature (K) and gas composition (mole %).
28
Page 28 of 38
P Pin and T Tin
Solve:
SRK-EOS:
an
G f (Pin , Tin )
Film mass continuity Film mom. continuity Gas mass continuity Gas mom. continuity
M
Solve:
us
Input:
cr
ip t
Figure_1
-
Eq. (2) Eq. (4) Eq. (1) Eq. (3)
ed
Calculate: G f ( momentum , mass transfer ) SRK-EOS:
pt
Solve:
Ac
ce
Solve:
new old Pnew f (G ,T )
Film Enthalpy - Eq. (6) Gas Enthalpy - Eq. (5)
ti t i 1 t yes
t i t END no END
Page 29 of 38
an
us
cr
ip t
Figure_2
Ac
ce pt
ed
M
Fig. 2
Page 30 of 38
M
an
us
cr
ip t
Figure_3
Ac
ce pt
ed
Fig. 3
Page 31 of 38
Figure_4
0 PT-2 Measured PT-2 Predicted PT-8 Measured PT-8 Predicted
-5
-15 -20
ip t
-25 -30 -35 10
20
30 Time (ms)
40
50
an
us
-40
cr
o
Temperature ( C )
-10
Ac
ce pt
ed
M
Fig. 4
Page 32 of 38
Figure_5
10
PT-8 Measured PT-8 Predicted PT-10 Measured PT-10 Predicted
8
ip t
7 6 5 40 60 Time (ms)
80
an
us
20
cr
Pressure ( MPa )
9
Ac
ce pt
ed
M
Fig. 5
Page 33 of 38
Figure_6
1.0
0.8 0.7 0.6
Measured 2 Pred. (MC= 12.25 kg/m /s)
0.5
Pred. (MC= 49 kg/m /s)
ip t
Pressure Ratio ( - )
0.9
2
2
0
50 100 150 200 250 300 Decompression Wave Speed (m/s)
350
an
us
0.4
cr
Pred. (MC= 6.125 kg/m /s)
Ac
ce pt
ed
M
Fig. 6
Page 34 of 38
Figure_7
3.5
Liquid film Gas (9.8 MPa) Gas (9.0 MPa) Gas (8.0 MPa) Gas (7.0 MPa) Gas (6.0 MPa)
2.5 2.0
ip t
1.5 1.0 0.5 0
2
4
6 Time (ms)
8
10
an
us
0.0
cr
Pipe Axial Elevation ( m )
3.0
Ac
ce pt
ed
M
Fig. 7
Page 35 of 38
Figure_8
80 Predicted - z = 0.1 m Predicted - z = 0.5 m Predicted - z = 1.0 m
70
50 40 30
ip t
Slip Velocity ( m/s )
60
20 10
2
4 6 Time (ms)
8
10
an
us
0
cr
0
Ac
ce pt
ed
M
Fig. 8
Page 36 of 38
Figure_9
120
Gas ( z = 0.1 m) Liquid ( z = 0.1 m) Gas ( z = 0.5 m) Liquid ( z = 0.5 m) Gas ( z = 1.0 m) Liquid ( z = 1.0 m)
100
60 40
ip t
Velocity (m/s)
80
0 1
2
3
4 5 6 Time (ms)
7
8
9
10
an
us
0
cr
20
Ac
ce pt
ed
M
Fig. 9
Page 37 of 38
Figure_10
11 z = 0.1 m z = 0.5 m z = 1.0 m
10
8 7 6
ip t
Pressure (MPa)
9
5
0
1
2
3
4 5 6 Time (ms)
7
8
9
10
an
us
3
cr
4
Ac
ce pt
ed
M
Fig. 10
Page 38 of 38