- Email: [email protected]
Mathematical modelling on rapid decompression in base natural gas mixtures under rupturing
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Contents lists available at SciVerse ScienceDirect
Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd
Mathematical modelling on rapid decompression in base natural gas mixtures under rupturing E. Burlutskiy ∗ A*STAR Institute of High Performance Computing, 1 Fusionopolis Way, 138632, Singapore
a b s t r a c t The paper presents a one-dimensional mathematical model of transient compressible thermal multi-component gas mixture flows in pipes. The set of the mass, momentum and enthalpy conservation equations is solved for gas phase in the model. Thermo-physical properties of multi-component gas mixture are calculated by solving the Equation of State (EOS) model. The Soave–Redlich–Kwong (SRK–EOS) model is chosen. Gas mixture viscosity is calculated on the basis of the Lee–Gonzales–Eakin (LGE) correlation. Numerical analysis on rapid decompression process in a shock tube having base natural gases is performed by using the proposed mathematical model. The model is successfully validated on the experimental measurements of the decompression wave speed in base natural gas mixtures. The proposed mathematical model shows a very good agreement with the experiments in a wide range of pressure values and predicts the decompression in base natural gases much better than analytical and mathematical models, which are available from the open source literature. © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. Keywords: Mathematical; Modelling; Rapid; Gas mixture; Decompression
1.
Introduction
A rupture on the pipeline wall due to fracture brings numerous problems for oil and gas engineers. It takes a lot of efforts fixing the problem. The fracture on the pipe wall does not propagate far and is arrested quickly, if a pipeline has been designed by employing of new advances in fracture propagation control. The fracture propagation control in natural gas pipeline transport service is usually made by using the Battelle two-curve method, which was developed by the Battelle Columbus Laboratories in order to determine the fracture arrest toughness (Eiber et al., 1993, 2004). The fracture propagation speed in the pipeline wall and the decompression wave speed in gas mixtures are used together in the Battelle analysis. The fracture propagation is arrested, if the decompression wave speed in gas mixture is quicker than the fracture propagation velocity in the pipeline wall material. Therefore, the information about the decompression wave speed in different gas mixtures is very important for the fracture propagation control and for the pipeline design. The fluid mixture composition, pipeline inner diameter, pressure, and temperature significantly influence on
∗
the decompression wave speed. The decompression in natural gas mixtures is very rapid non-isothermal process. A transient mathematical model of compressible thermal multi-component gas mixture flow in pipes gives very detailed information on basic flow parameter’s distribution and flow behaviour in a wide range of operating parameters. Information about the mathematical modelling on rapid decompression process in natural gas mixtures is extremely limited in the open source literature. Most part of papers contains very simplified engineering analysis. Experimental measurements of the decompression wave speed in dry, base and rich natural gas mixtures were conducted by TransCanada PipeLines (TCPL) (Botros et al., 2003, 2004, 2007, 2010a,b). Measurements were performed much more intensively compared to mathematical modelling. The influence of a shock tube inner diameter, gas mixture composition, pressure, and temperature was examined in details experimentally. Pressure values in those studies were varied in the range between 10 MPa and 37 MPa (Botros et al., 2003, 2010a,b). Most of measurements were made on small-diameter shock tubes, where the friction force influences on the flow behaviour much stronger compared to large-diameter pipes.
Tel.: +65 6419 1386; fax: +65 6467 4350. E-mail addresses: [email protected], [email protected] Received 15 December 2011; Received in revised form 8 May 2012; Accepted 25 June 2012 0263-8762/$ – see front matter © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cherd.2012.06.017
64
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Nomenclature Dpipe diameter of the pipe (m) D1 , D2 , D3 parameters of the viscosity correlation DP1 , DP8 , DP5 , DP6 distances between PT locations and the rupture end of the pipe (m) enthalpy of the fluid (J/kg) hG molar mass of the gas mixture (g/mol) M molecular weight of the gas mixture (g/mol) MWg N number of components in the gas mixture P total pressure (MPa) critical pressure value of the component i of gas PCi mixture (MPa) Pinitial initial pressure in the shock tube before rupturing (MPa) R universal gas constant (J/mol/K) RG−Wall gas mixture to wall friction term S cross-sectional area of the pipe (m2 ) time (s) t T temperature of the gas mixture (K) TCi critical temperature value of the component i of gas mixture (K) Tinitial initial temperature in the shock tube before rupturing (K) UG velocity of the gas mixture (m/s) V volume of the gas mixture (m3 ) axial co-ordinate of the decompression tube z (m) zi mole fraction of the component i of gas mixture Z compressibility factor of the gas mixture s volume fraction of the gas mixture internal surface roughness (m) ε gas mixture viscosity (Pa s) G G−Wall gas to wall friction coefficient ˘ perimeter of the pipe (m) G density of the gas mixture (kg/m3 ) G−Wall gas mixture to wall friction term acentric factor of the component i of gas mixωi ture
The program called GASDECOM (Eiber et al., 1993) is the most well-known software in the field of natural gas transport engineering. It calculates the decompression wave speed values in natural gases (Botros et al., 2003, 2007) by using analytical correlations for gas mixture hydrodynamics together with different EOS models. The friction force is not taken into account in GASDECOM’s model. Thus, the application area is limited on large-diameter pipelines, where accurate calculation of thermo-physical fluid properties is much more important compared to proper modelling of fluid hydrodynamic behaviour. The program predicts the decompression wave speed values with a reasonably good level of accuracy compared to the experimental data if those experimental values are determined from pressure transducers, which are mounted in the rupture end area of the tube, where the friction’s influence is not very strong. The comparison between measured data and GASDECOM calculations is poor, if the decompression wave speed is determined experimentally from pressure transducers, which are located far away from the rupture end of the pipe, and where the friction force influences on the flow behaviour significantly. Numerical simulations on rapid decompression process in rich and base gas
mixtures were performed (Botros et al., 2007) by using the commercial one-dimensional OLGA code (SPT-group) (Bendiksen et al., 1991) as well. All predictions, which were made by using OLGA code (Botros et al., 2007), show a poor comparison with experimental data and all calculated values are significantly over-predicted. The paper presents a one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flow in pipes. Numerical analysis on rapid decompression process in a shock tube having base natural gases is performed by using the proposed mathematical model. The model is successfully validated on the experimental data on rapid decompression in base natural gas mixtures (Botros et al., 2003) and it shows a good agreement with the experiments. A proper modelling of gas mixture continuity, momentum and enthalpy equations together with taking into account the thermo-physical fluid properties (SRK–EOS) is the key to successful model performance. The proposed mathematical model predicts the decompression wave speed in base natural gases much better compared to mathematical and analytical models, which are available from the open source literature.
2. One-dimensional mathematical model of transient thermal fluid flow in pipes The set of mass, momentum and enthalpy conservation equations is solved for gas phase in the mathematical model. This set of equations for single phase gas mixture flow in general form is written as (Wallis, 1969): ∂˛G G UG ∂˛G G + =0 ∂t ∂z
(1)
2 ∂˛G G UG ∂˛G G UG ∂P + = −˛G − RG−Wall ∂t ∂z ∂z
(2)
∂˛G G UG hG ∂˛G G hG + = ˛G ∂t ∂z
∂P ∂t
+ UG
∂P ∂z
+ RG−Wall UG
(3)
Here, ˛G is the volume fraction of the gas mixture; G is the density of the gas mixture; UG is the velocity of the gas mixture; P is the total pressure; RG−Wall is the friction term, hG is the enthalpy of the fluid, t is the time, z is the axial co-ordinate. The last term in the equation of enthalpy conservation (3) represents an additional contribution into the system’s enthalpy due to the friction force. The friction term is written in the form of Blasius (1913):
RG−Wall =
˘ , S G−Wall
G−Wall = 64/ReG ,
G−Wall =
G UG Dpipe G
8
ReG < 1600
G−Wall = 0.316/Re0.25 G , ReG =
2 G−Wall G UG
ReG > 1600
(4)
(5)
(6)
Here, ˘ is the perimeter of the pipe; S is the cross-sectional area of the pipe; G−Wall is the friction term, which represents friction between the gas and the wall of the pipe; G−Wall is the friction coefficient; Dpipe is the diameter of the pipe; G is the viscosity of the fluid.
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
3. Thermo-physical properties of natural gas mixture Thermo-physical fluid properties are modelled by solving of the Equation of State (EOS) in the form of the Soave–Redlich–Kwong (Soave, 1979) model. The set of equations and correlations (SRK–EOS) may be written as (Soave, 1979): P=
a=
RT a − (V − b) V(V + b) N N
zi zj
(7)
ai aj (1 − kij );
b=
i=1 j=1
ai = 0.42748
N
zi bi
(8)
i=1
2 R2 TCi
PCi
1 + mi
1−
mi = 0.48 + 1.574 ωi − 0.176ωi2
T TCi
2 ;
bi = 0.08664
RTCi PCi (9)
(10)
Here, V is the volume of the gas mixture; N is the number of components in the gas mixture; T is the temperature of the gas mixture; R is the universal gas constant; ωi is the acentric factor of the component i; PCi , TCi are critical values of the pressure and temperature, correspondently; zi is the mole fraction of the component i. The compressibility factor (Z) of the gas mixture is calculated from the following equation (Soave, 1979): Z3 − Z2 + (A − B − B2 )Z − AB = 0 A=
aP ; R2 T 2
P = G
B=
bP RT
(11) (12)
R TZ M
(13)
Here, M is the molar mass of a substance. The viscosity of the gas mixture is calculated by using of the Lee–Gonzales–Eakin (LGE) correlation and may be written as (Lee et al., 1966): D3 ) G = D1 × 10−4 · exp (D2 G
(14)
(9.38 + 0.016 MWg )T 1.5 209.2 + 19.26 MWg + T
(15)
D1 =
D2 = 3.448 +
986.4 T
D3 = 2.447 − 0.224 · D2
+ 0.01 MWg
(16) (17)
Here, MWg is the molecular weight of the gas mixture; g is the gas density in (g/cc) here; T is the temperature in (R) here. The algorithm of solving 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 known as the Thomas algorithm (Patankar, 1980). It is a simplified form of Gaussian elimination that can be used to solve tri-diagonal systems of equations. The set of unsteady governing equations is transformed into the standard form of the discrete
65
analog of the tri-diagonal system (Patankar, 1980) by using the fully implicit numerical scheme (Burlutskiy, 2012). In this case the equation is reduces to the steady state discretization equation, if the time step goes to infinity. The following numerical scheme is applied. The computational pipe is divided into N parts, which are called the control volumes. The set of governing equations is transformed into the discrete analog of the tri-diagonal system, which is calculated at each control volume of the computational pipe. The velocity of gas is set-up to zero (Ui = 0 m/s; i = 1, . . ., N; N is the number of control volumes of the pipe) in each control volume of the computational pipe once at the time before rupture started. The pressure and temperature of gas mixture is setup to initial values before rupturing (Ti = Tinitial ; Pi = Pinitial ; i = 1, . . ., N) in each control volume of the computational pipe. The pressure outside the pipe (i.e. after the rupture end of the tube) is set-up to the atmospheric value (PN+1 = 0.1 MPa). Thus, the pressure gradient is non-zero at the rupture end of the pipe (PN = Pinitial ; PN+1 = 0.1 MPa; PN PN+1 ). The density of gas mixture is calculated first time (before rupturing) by using the Soave–Redlich–Kwong equation of State (SRK–EOS) model. This is the direct way of the SRK–EOS use. The density is a function of the pressure and temperature. It is written as: initial = f (Pinitial , Tinitial )
(18)
The density of gas mixture is set-up then to calculated values (i = initial ; i = 1, . . ., N) in each control volume of the computational pipe. The volume fraction of gas mixture is set-up to 1 (one) in each control volume of the computational pipe at the beginning. The volume fraction is conceived as a local variable due to compressible nature of the problem. The value of this variable is changed during the time step from 1 (one) at the beginning (before each time step) to other value (i.e. between one and zero) after the time step. This change is found by solving the continuity conservation equation, where the density is a constant (during each time step only) and the volume fraction is a variable during the time step. When the time step is made and a change in volume fraction is calculated, the density is determined by using new volume fraction values, and the volume fraction is set-up to 1 (one) at each control volume of the pipe again. Thus, the volume fraction of gas is set-up to 1 (one) at each control volume of the pipe before every time step. The momentum conservation equation is solved for the first time step. The velocity of gas mixture is calculated in each control volume of the pipe. The pressure gradient at the rupture end of the pipe is a driving force for the momentum in this case. The continuity conservation equation is solved after the momentum equation. The density of gas is assumed to be constant during each time step in each control volume. Therefore, volume fraction values are calculated for each control volume of the pipe by solving the continuity equation, which is influenced by the momentum. Local iterations between continuity and momentum equations are performed at each time step in order to get a better convergence and a stability of the numerical scheme. New values of gas density are calculated by using new values of gas volume fraction, which are found by solving both momentum and continuity equations. Since new density values are obtained, the volume fraction values of gas mixture are set-up to 1 (one) at each control volume of the computational pipe again. A continuous change in gas
66
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Fig. 1 – Schematic of the experimental decompression tube. density is provided by small value of the time step. New pressure values are calculated thereafter by using the SRK–EOS model determining the inverse correspondence as: P = f (G , T)
(19)
New pressure values are obtained then at each control volume of the pipe. The pressure value outside the pipe (i.e. after the rupture end of the tube) is constantly set-up to the atmospheric value (PN+1 = 0.1 MPa) at each time step. Temperature values are calculated thereafter by solving the enthalpy conservation equation in each control volume of the computational pipe. Thus, the density, pressure, temperature and velocity values are new (i.e. re-calculated) now. Then this procedure is repeated again for the new time step. A one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flow in pipes was developed (Burlutskiy, 2012) under the research project named “Multi-Component Natural Gas mixture flows in Pipelines and Wells” in PETROSOFT-DC. This mathematical model was implemented into the FORTRAN computer code and was named the “Gas Decompression Program” (GDP code).
4. Numerical analysis on rapid decompression process in base natural gas mixtures The proposed mathematical model is validated on the experimental data on rapid decompression in a shock tube having base natural gas mixtures (Botros et al., 2003). Rich gas mixtures are not considered in this study, because a liquid condensation from gaseous phase is usually occurred at rich and medium rich gases. The presented model does not account for a liquid condensation phenomenon in this study. The model will be extended into the case of two-phase gas–liquid flows in pipes, where a liquid condensation will be taken into account, and the model will be validated on the experimental measurements of rapid decompression in pipe systems having rich gas mixtures (Botros et al., 2007) in other work. Results of this work will be presented in future papers.
The experimental measurements of the decompression wave speed in natural gases were conducted by TCPL (TransCanada PipeLines) at TCPL Gas Dynamic Test Facility in Didsbury, Alberta, Canada (Botros et al., 2003). The main test section of the facility is a shock tube having a length of 30 m and the inner diameter of 49.325 mm. The internal surface of the shock tube has a roughness better than 40 micro-inches (i.e. less than 1.016 × 10−6 m). Basing on the value of the roughness-to-diameter ratio (ε/Dpipe = 20.5 × 10−6 ; ε = 1.016 × 10−6 m; Dpipe = 49.325 × 10−3 m), the internal surface of the shock tube is determined as the “smooth” tube surface (Botros et al., 2010b). Therefore, the wall roughness does not influence on the flow parameters strongly. The classical Blasius correlation (4)–(6) for the smooth pipes is allowed in this case. The influence of the internal surface roughness will be examined numerically in details by using different friction correlations, where the internal roughness will be taken into account, in other future work as well as other experimental measurements (Botros et al., 2010b) will be employed for the model validation. A rupture disc is placed at one end of the shock tube, which is upon rupturing. A decompression wave propagates up into the pressurized test section. A few high frequency responses pressure transducers (PT) are mounted into the tube in order to capture the time history values of the expansion fan (Botros et al., 2003). Several rupture discs were introduced into the shock tube end for different pressures at rupture. Decompression wave speed values are determined from the time between signals from P1 to P8 pressure transducers as well as the time between signals from P5 to P6 transducers. Fig. 1 shows the schematic view of the experimental decompression tube. Despite the fact that the main test section of the shock tube has a length of 30 m, the computational decompression pipe having a length of 50 m is simulated from the point of view of the numerical scheme stability. The inner diameter of the computational pipe is 49.325 mm. One end of the shock tube is selected to be the closed end. The rupture disc boundary condition is modelled at other pipe end. Distances between the rupture disc and pressure transducers are selected to be identical to real distances on the experimental
Table 1 – Initial pressure (MPa), temperature (K) and gas composition (mole %). Cases
Pinitial
Tinitial
N2
CO2
C1
C2
C3
i-C4
n-C4
i-C5
n-C5
Case 1 Case 2 Case 3
10.41 13.8 20.67
274.07 264.77 264.72
0.697 0.647 0.699
1.097 1.197 1.279
92.955 93.02 92.757
4.076 3.876 4.075
0.8 0.909 0.861
0.099 0.108 0.103
0.137 0.158 0.146
0.066 0.059 0.053
0.073 0.026 0.027
67
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Pressure (MPa)
9 8 7 6 5
14
PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
12
Pressure (MPa)
PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
10
10
8
6
4 0
2
4
6
8 10
40
50
60
70
80
90
100
Time (ms)
decompression tube as it is shown in Fig. 1. Those distances between PT P1, P8, P5, P6 and the rupture disc are 0.4, 1.15, 16.4, 18.4 m, correspondently. Three experimental measurement data cases (Botros et al., 2003) having different initial values are selected in order to be simulated (Table 1). The following gas mixture composition (base natural gases), initial pressure and initial temperature are selected (Botros et al., 2003). Predictions are started with the initial pressure of 10.58 MPa and the temperature of 247.4 K in each computational cell of the pipe for the case 1 (Table 1). New values of the velocity, temperature, density and pressure are re-calculated after each time step. The upper limit of the time step is selected from the point of view of the numerical stability of calculations. The decompression is very quick process and basic parameters of gas mixture are changed very rapidly in time. Therefore, the time step has to be small enough in order to simulate continues change in pressure and temperature values in time at every computational cell. The computational time step is equal to 2 × 10−5 s in all of considered cases. Pressure values are collected at P1, P8, P5 and P6 PT locations after each time step. Fig. 2 shows the pressure time history values at P1, P8, P5, P6 PT locations for the case 1. Experimental points are shown in all figures as symbols. Continues lines represent predictions, which are made by using the proposed model. Fig. 3 shows measured (Botros et al., 2003)
1.0 Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
Pressure ratio (-)
0.8 0.7 0.6 0.5 50
100
150
200
250
5
10
15
40
60
80
100
120
140
Time (ms)
Fig. 2 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 1.
0.9
0
300
350
400
Decompression Wave Speed (m/s) Fig. 3 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 1.
Fig. 4 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 2.
and predicted values of the decompression wave speed, which are determined from P1 to P8 and P5 to P6 PT locations, correspondently. Pressure values are normalized on the initial pressure value before rupturing. Figs. 2–3 show a good agreement between experimental data and predictions, which are made by using GDP code. The proposed model simulates the slope in the pressure time history values correctly at each PT location and shows a good agreement with the measured data. The pressure shows rapid decrease at P1 and P8 PT locations. Those PT are located near the rupture end of the pipe. The friction does not influence on the decompression wave speed strongly here and, therefore, it does not slow it down significantly. The slope of the pressure time history curves at P5 and P6 PT locations is more flatter compared to P1 and P8 due to strong influence of the gas to wall friction on gas mixture release from the pipe. Calculations of the decompression wave speed in base natural gas mixtures in considered experimental conditions (Table 1; Cases 1–3) were made by TCPL (Botros et al., 2003) by using the analytical GASDECOM (Eiber et al., 1993) code. In order to compare those predictions with GDP calculations, GASDECOM’s calculated values are taken from Botros et al. (2003) and are put into Figs. 3, 5 and 7 of the paper. Broken curves in those figures represent GASDECOM’s numerical data (Botros et al., 2003). GASDECOM’s predictions are in good agreement with the proposed model and experimental decompression wave speed values, which are determined from P1 and P8 PT locations. However, GASDECOM does not able to predict the decompression wave speed values, which are determined from P5 and P6 PT locations, because the friction force is not taken into account in GASDECOM’s analysis. The proposed model (i.e. GDP code) shows ability in calculating the decompression wave speed, which is determined from P5 to P6 PT locations, correctly and with high order of agreement with the measurements. Other predictions of the decompression wave speed values are made by using the same shock tube topology (Fig. 1), but having much higher initial pressure before rupture started (Table 1, Case 2). The initial temperature and gas composition are up-dated accordingly to the experimental values, which were selected for the simulations (Table 1). Pressure time history values show a good agreement with experimental data in this case (Figs. 4–5), also. Predictions, which are made by using GDP code, show much better agreement
68
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
1.0
1.0 Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
0.8
Pressure ratio (-)
Pressure ratio (-)
0.9
0.9
0.7 0.6
0.7 0.6 0.5 0.4
0.5 0.4
0.8
Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
50
100
150
200
250
300
350
0.3 100
400
200
300
400
500
Decompression Wave Speed (m/s)
Decompression Wave Speed (m/s) Fig. 5 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 2.
Fig. 7 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 3.
with the experimental data (Fig. 5) compare to GASDECOM’s calculations, which are taken from Botros et al. (2003) again. The last set of predictions of the decompression wave speed are made for the same shock tube topology (Fig. 1) but having the highest initial pressure value (Table 1, Case 3). The agreement between predicted pressure time history values and experimental measurements has a high order of magnitude (Fig. 6). The proposed model simulates the decompression process in base natural gas mixtures correctly and provides very detailed information about local distribution of basic flow parameters (spatial and time history, both), which is very helpful in engineering analysis. The prediction, which is made by using the proposed model, shows much better agreement with the experimental data compare to analytical GASDECOM model (Fig. 7). The analytical GASDECOM (Eiber et al., 1993) is one of the most frequently used engineering software in gas field applications performing quick simulations of the decompression wave speed in natural gases. However, this model does not account for the friction between the gas and the pipe’s wall. The flow behaviour in a pipe’s area far away from the rupture place is totally different compared to the rupture end area due to friction force contribution. It was observed experimentally, that the decompression wave speed, which is determined from P5 to P6 PT locations, farther away from the rupture end,
is much lower than the corresponding one (Figs. 2–7), which is determined from a closer pair of P1 to P8 PT locations (Botros et al., 2003). Therefore, the analytical GASDECOM program itself does not simulate the decompression process in natural gas mixtures correctly in the case of small-diameter shock tubes, where the friction between the gas and the wall plays an important role. GASDECOM’s values are over-predicted in this case. Its application area is limited on the large-diameter pipeline cases due to this fact. Simulations on rapid decompression process in dry and base gas mixtures were made by TCPL (Botros et al., 2007) by using the well-know and frequently used commercial 1D software OLGA, which is developed by SPT-group (Bendiksen et al., 1991). All predictions, which were made by using OLGA, are found in bad agreement with the experimental data. The frontal wave speed values, which are predicted by OLGA code, are significantly lower than the measured one (Botros et al., 2007). Thus, the proposed mathematical model of transient compressible thermal multi-component gas mixture flows in pipes predicts the decompression process in a shock tube having base natural gases much better than analytical and mathematical models, which are available from the open source literature. Moreover, simulations, which are made by using the presented mathematical model, are quick in time and allow simulating of large-scale industrial pipes and shock tubes.
22 PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
20
Pressure (MPa)
18 16 14 12 10 8 6 0 2 4 6 8 10
30 40
50 60 70
80 90
Time (ms) Fig. 6 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 3.
5.
Conclusions
A one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flows in pipes is presented in the paper. The set of mass, momentum and enthalpy conservation equations for gas phase is solved in the model. Thermo-physical properties of multi-component gas mixture are calculated by solving the Equation of State (EOS) model. The Soave–Redlich–Kwong (SRK–EOS) model is chosen. Gas mixture viscosity is calculated on the basis of the Lee–Gonzales–Eakin (LGE) correlation. Numerical analysis on rapid decompression process in base natural gas mixtures, which is performed by using the proposed mathematical model, is presented in the paper. The model is successfully validated on the experimental data on rapid decompression in a shock tube having base natural gas mixtures from. Three experimental cases with high, medium
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
and low initial pressure values before rupturing are simulated. The proposed mathematical model shows a very good agreement with the experiments in a wide range of pressure values. It has been shown, that the mathematical model predicts the decompression in base natural gases much better than analytical and mathematical models, which are available from the open source literature. The presented model is highly necessary and useful in pipeline design and flow assurance. The minimum of fracture arrest toughness of the pipe wall material may be determined on the basis of the Battelle two-curve method with taking into account of the proposed model together with fracture propagation speed model. The influence of the pressure, temperature, fluid composition, and pipeline diameter may be examined by using the presented model.
Acknowledgement The author would like to acknowledge the support of PETROSOFT-DC for giving possibility to use the computer program (GDP code) performing the mathematical modelling and numerical analysis on rapid decompression in a shock tube having base natural gases.
References Bendiksen, K.H., Maines, D., Moe, R., Nuland, S., 1991. The dynamic two-fluid model OLGA: theory and application. SPE Prod. Eng. V6 (N2), 171–180. Blasius, P.R.H., 1913. Das Aehnlichkeitsgesetz bei Reibungsvorgangen in Fluessigkeiten. 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. In: American Gas Association Proceedings (AGA-2003).
69
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2004. Determination of decompression wave speed in rich gas mixtures. Can. J. Chem. Eng. 82, 880–891. 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. Pres. Ves. Pip. 84, 358–367. 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. Pres. Ves. Technol. 132 (051303), 1–15. Botros, K.K., Geerligs, J., Rothwell, B., Carlson, L., Fletcher, L., Venton, P., 2010b. 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. Pres. Ves. Pip. 87, 681–695. Burlutskiy, E., 2012. Mathematical modelling of non-isothermal multi-component fluid flow in pipes applying to rapid gas decompression in rich and base natural gases. In: Proceedings of the International Conference on Fluid Mechanics, Heat Transfer and Thermodynamics–ICFMHTT-2012, vol. 61, Zurich, Switzerland, 15–17 January, pp. 156–167. Eiber, R.J., Bubenik, T.A., Maxey, W.A., 1993. GASDECOM, Computer Code for the Calculation of Gas Decompression Speed. Fracture Control for Natural Gas Pipelines. PRCI Report No. L51691. Eiber, R.J., Carlson, L., Leis, B., 2004. Fracture control requirements for gas transmission pipelines. In: Proceedings of the Fourth International Conference on Pipeline Technology, p. 437. Lee, A.L., Gonzalez, M.N., Eakin, B.E., 1966. The viscosity of natural gases. J. Petrol. Technol., 997–1000. Patankar, S., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York. Soave, G., 1979. Equilibrium constants from a modified Redlich–Kwong equation of state. Chem. Eng. Sci. 27, 1197–1203. Wallis, G.B., 1969. One-dimensional Two-phase Flows. McGraw Hill, New York.
Contents lists available at SciVerse ScienceDirect
Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd
Mathematical modelling on rapid decompression in base natural gas mixtures under rupturing E. Burlutskiy ∗ A*STAR Institute of High Performance Computing, 1 Fusionopolis Way, 138632, Singapore
a b s t r a c t The paper presents a one-dimensional mathematical model of transient compressible thermal multi-component gas mixture flows in pipes. The set of the mass, momentum and enthalpy conservation equations is solved for gas phase in the model. Thermo-physical properties of multi-component gas mixture are calculated by solving the Equation of State (EOS) model. The Soave–Redlich–Kwong (SRK–EOS) model is chosen. Gas mixture viscosity is calculated on the basis of the Lee–Gonzales–Eakin (LGE) correlation. Numerical analysis on rapid decompression process in a shock tube having base natural gases is performed by using the proposed mathematical model. The model is successfully validated on the experimental measurements of the decompression wave speed in base natural gas mixtures. The proposed mathematical model shows a very good agreement with the experiments in a wide range of pressure values and predicts the decompression in base natural gases much better than analytical and mathematical models, which are available from the open source literature. © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. Keywords: Mathematical; Modelling; Rapid; Gas mixture; Decompression
1.
Introduction
A rupture on the pipeline wall due to fracture brings numerous problems for oil and gas engineers. It takes a lot of efforts fixing the problem. The fracture on the pipe wall does not propagate far and is arrested quickly, if a pipeline has been designed by employing of new advances in fracture propagation control. The fracture propagation control in natural gas pipeline transport service is usually made by using the Battelle two-curve method, which was developed by the Battelle Columbus Laboratories in order to determine the fracture arrest toughness (Eiber et al., 1993, 2004). The fracture propagation speed in the pipeline wall and the decompression wave speed in gas mixtures are used together in the Battelle analysis. The fracture propagation is arrested, if the decompression wave speed in gas mixture is quicker than the fracture propagation velocity in the pipeline wall material. Therefore, the information about the decompression wave speed in different gas mixtures is very important for the fracture propagation control and for the pipeline design. The fluid mixture composition, pipeline inner diameter, pressure, and temperature significantly influence on
∗
the decompression wave speed. The decompression in natural gas mixtures is very rapid non-isothermal process. A transient mathematical model of compressible thermal multi-component gas mixture flow in pipes gives very detailed information on basic flow parameter’s distribution and flow behaviour in a wide range of operating parameters. Information about the mathematical modelling on rapid decompression process in natural gas mixtures is extremely limited in the open source literature. Most part of papers contains very simplified engineering analysis. Experimental measurements of the decompression wave speed in dry, base and rich natural gas mixtures were conducted by TransCanada PipeLines (TCPL) (Botros et al., 2003, 2004, 2007, 2010a,b). Measurements were performed much more intensively compared to mathematical modelling. The influence of a shock tube inner diameter, gas mixture composition, pressure, and temperature was examined in details experimentally. Pressure values in those studies were varied in the range between 10 MPa and 37 MPa (Botros et al., 2003, 2010a,b). Most of measurements were made on small-diameter shock tubes, where the friction force influences on the flow behaviour much stronger compared to large-diameter pipes.
Tel.: +65 6419 1386; fax: +65 6467 4350. E-mail addresses: [email protected], [email protected] Received 15 December 2011; Received in revised form 8 May 2012; Accepted 25 June 2012 0263-8762/$ – see front matter © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cherd.2012.06.017
64
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Nomenclature Dpipe diameter of the pipe (m) D1 , D2 , D3 parameters of the viscosity correlation DP1 , DP8 , DP5 , DP6 distances between PT locations and the rupture end of the pipe (m) enthalpy of the fluid (J/kg) hG molar mass of the gas mixture (g/mol) M molecular weight of the gas mixture (g/mol) MWg N number of components in the gas mixture P total pressure (MPa) critical pressure value of the component i of gas PCi mixture (MPa) Pinitial initial pressure in the shock tube before rupturing (MPa) R universal gas constant (J/mol/K) RG−Wall gas mixture to wall friction term S cross-sectional area of the pipe (m2 ) time (s) t T temperature of the gas mixture (K) TCi critical temperature value of the component i of gas mixture (K) Tinitial initial temperature in the shock tube before rupturing (K) UG velocity of the gas mixture (m/s) V volume of the gas mixture (m3 ) axial co-ordinate of the decompression tube z (m) zi mole fraction of the component i of gas mixture Z compressibility factor of the gas mixture s volume fraction of the gas mixture internal surface roughness (m) ε gas mixture viscosity (Pa s) G G−Wall gas to wall friction coefficient ˘ perimeter of the pipe (m) G density of the gas mixture (kg/m3 ) G−Wall gas mixture to wall friction term acentric factor of the component i of gas mixωi ture
The program called GASDECOM (Eiber et al., 1993) is the most well-known software in the field of natural gas transport engineering. It calculates the decompression wave speed values in natural gases (Botros et al., 2003, 2007) by using analytical correlations for gas mixture hydrodynamics together with different EOS models. The friction force is not taken into account in GASDECOM’s model. Thus, the application area is limited on large-diameter pipelines, where accurate calculation of thermo-physical fluid properties is much more important compared to proper modelling of fluid hydrodynamic behaviour. The program predicts the decompression wave speed values with a reasonably good level of accuracy compared to the experimental data if those experimental values are determined from pressure transducers, which are mounted in the rupture end area of the tube, where the friction’s influence is not very strong. The comparison between measured data and GASDECOM calculations is poor, if the decompression wave speed is determined experimentally from pressure transducers, which are located far away from the rupture end of the pipe, and where the friction force influences on the flow behaviour significantly. Numerical simulations on rapid decompression process in rich and base gas
mixtures were performed (Botros et al., 2007) by using the commercial one-dimensional OLGA code (SPT-group) (Bendiksen et al., 1991) as well. All predictions, which were made by using OLGA code (Botros et al., 2007), show a poor comparison with experimental data and all calculated values are significantly over-predicted. The paper presents a one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flow in pipes. Numerical analysis on rapid decompression process in a shock tube having base natural gases is performed by using the proposed mathematical model. The model is successfully validated on the experimental data on rapid decompression in base natural gas mixtures (Botros et al., 2003) and it shows a good agreement with the experiments. A proper modelling of gas mixture continuity, momentum and enthalpy equations together with taking into account the thermo-physical fluid properties (SRK–EOS) is the key to successful model performance. The proposed mathematical model predicts the decompression wave speed in base natural gases much better compared to mathematical and analytical models, which are available from the open source literature.
2. One-dimensional mathematical model of transient thermal fluid flow in pipes The set of mass, momentum and enthalpy conservation equations is solved for gas phase in the mathematical model. This set of equations for single phase gas mixture flow in general form is written as (Wallis, 1969): ∂˛G G UG ∂˛G G + =0 ∂t ∂z
(1)
2 ∂˛G G UG ∂˛G G UG ∂P + = −˛G − RG−Wall ∂t ∂z ∂z
(2)
∂˛G G UG hG ∂˛G G hG + = ˛G ∂t ∂z
∂P ∂t
+ UG
∂P ∂z
+ RG−Wall UG
(3)
Here, ˛G is the volume fraction of the gas mixture; G is the density of the gas mixture; UG is the velocity of the gas mixture; P is the total pressure; RG−Wall is the friction term, hG is the enthalpy of the fluid, t is the time, z is the axial co-ordinate. The last term in the equation of enthalpy conservation (3) represents an additional contribution into the system’s enthalpy due to the friction force. The friction term is written in the form of Blasius (1913):
RG−Wall =
˘ , S G−Wall
G−Wall = 64/ReG ,
G−Wall =
G UG Dpipe G
8
ReG < 1600
G−Wall = 0.316/Re0.25 G , ReG =
2 G−Wall G UG
ReG > 1600
(4)
(5)
(6)
Here, ˘ is the perimeter of the pipe; S is the cross-sectional area of the pipe; G−Wall is the friction term, which represents friction between the gas and the wall of the pipe; G−Wall is the friction coefficient; Dpipe is the diameter of the pipe; G is the viscosity of the fluid.
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
3. Thermo-physical properties of natural gas mixture Thermo-physical fluid properties are modelled by solving of the Equation of State (EOS) in the form of the Soave–Redlich–Kwong (Soave, 1979) model. The set of equations and correlations (SRK–EOS) may be written as (Soave, 1979): P=
a=
RT a − (V − b) V(V + b) N N
zi zj
(7)
ai aj (1 − kij );
b=
i=1 j=1
ai = 0.42748
N
zi bi
(8)
i=1
2 R2 TCi
PCi
1 + mi
1−
mi = 0.48 + 1.574 ωi − 0.176ωi2
T TCi
2 ;
bi = 0.08664
RTCi PCi (9)
(10)
Here, V is the volume of the gas mixture; N is the number of components in the gas mixture; T is the temperature of the gas mixture; R is the universal gas constant; ωi is the acentric factor of the component i; PCi , TCi are critical values of the pressure and temperature, correspondently; zi is the mole fraction of the component i. The compressibility factor (Z) of the gas mixture is calculated from the following equation (Soave, 1979): Z3 − Z2 + (A − B − B2 )Z − AB = 0 A=
aP ; R2 T 2
P = G
B=
bP RT
(11) (12)
R TZ M
(13)
Here, M is the molar mass of a substance. The viscosity of the gas mixture is calculated by using of the Lee–Gonzales–Eakin (LGE) correlation and may be written as (Lee et al., 1966): D3 ) G = D1 × 10−4 · exp (D2 G
(14)
(9.38 + 0.016 MWg )T 1.5 209.2 + 19.26 MWg + T
(15)
D1 =
D2 = 3.448 +
986.4 T
D3 = 2.447 − 0.224 · D2
+ 0.01 MWg
(16) (17)
Here, MWg is the molecular weight of the gas mixture; g is the gas density in (g/cc) here; T is the temperature in (R) here. The algorithm of solving 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 known as the Thomas algorithm (Patankar, 1980). It is a simplified form of Gaussian elimination that can be used to solve tri-diagonal systems of equations. The set of unsteady governing equations is transformed into the standard form of the discrete
65
analog of the tri-diagonal system (Patankar, 1980) by using the fully implicit numerical scheme (Burlutskiy, 2012). In this case the equation is reduces to the steady state discretization equation, if the time step goes to infinity. The following numerical scheme is applied. The computational pipe is divided into N parts, which are called the control volumes. The set of governing equations is transformed into the discrete analog of the tri-diagonal system, which is calculated at each control volume of the computational pipe. The velocity of gas is set-up to zero (Ui = 0 m/s; i = 1, . . ., N; N is the number of control volumes of the pipe) in each control volume of the computational pipe once at the time before rupture started. The pressure and temperature of gas mixture is setup to initial values before rupturing (Ti = Tinitial ; Pi = Pinitial ; i = 1, . . ., N) in each control volume of the computational pipe. The pressure outside the pipe (i.e. after the rupture end of the tube) is set-up to the atmospheric value (PN+1 = 0.1 MPa). Thus, the pressure gradient is non-zero at the rupture end of the pipe (PN = Pinitial ; PN+1 = 0.1 MPa; PN PN+1 ). The density of gas mixture is calculated first time (before rupturing) by using the Soave–Redlich–Kwong equation of State (SRK–EOS) model. This is the direct way of the SRK–EOS use. The density is a function of the pressure and temperature. It is written as: initial = f (Pinitial , Tinitial )
(18)
The density of gas mixture is set-up then to calculated values (i = initial ; i = 1, . . ., N) in each control volume of the computational pipe. The volume fraction of gas mixture is set-up to 1 (one) in each control volume of the computational pipe at the beginning. The volume fraction is conceived as a local variable due to compressible nature of the problem. The value of this variable is changed during the time step from 1 (one) at the beginning (before each time step) to other value (i.e. between one and zero) after the time step. This change is found by solving the continuity conservation equation, where the density is a constant (during each time step only) and the volume fraction is a variable during the time step. When the time step is made and a change in volume fraction is calculated, the density is determined by using new volume fraction values, and the volume fraction is set-up to 1 (one) at each control volume of the pipe again. Thus, the volume fraction of gas is set-up to 1 (one) at each control volume of the pipe before every time step. The momentum conservation equation is solved for the first time step. The velocity of gas mixture is calculated in each control volume of the pipe. The pressure gradient at the rupture end of the pipe is a driving force for the momentum in this case. The continuity conservation equation is solved after the momentum equation. The density of gas is assumed to be constant during each time step in each control volume. Therefore, volume fraction values are calculated for each control volume of the pipe by solving the continuity equation, which is influenced by the momentum. Local iterations between continuity and momentum equations are performed at each time step in order to get a better convergence and a stability of the numerical scheme. New values of gas density are calculated by using new values of gas volume fraction, which are found by solving both momentum and continuity equations. Since new density values are obtained, the volume fraction values of gas mixture are set-up to 1 (one) at each control volume of the computational pipe again. A continuous change in gas
66
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Fig. 1 – Schematic of the experimental decompression tube. density is provided by small value of the time step. New pressure values are calculated thereafter by using the SRK–EOS model determining the inverse correspondence as: P = f (G , T)
(19)
New pressure values are obtained then at each control volume of the pipe. The pressure value outside the pipe (i.e. after the rupture end of the tube) is constantly set-up to the atmospheric value (PN+1 = 0.1 MPa) at each time step. Temperature values are calculated thereafter by solving the enthalpy conservation equation in each control volume of the computational pipe. Thus, the density, pressure, temperature and velocity values are new (i.e. re-calculated) now. Then this procedure is repeated again for the new time step. A one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flow in pipes was developed (Burlutskiy, 2012) under the research project named “Multi-Component Natural Gas mixture flows in Pipelines and Wells” in PETROSOFT-DC. This mathematical model was implemented into the FORTRAN computer code and was named the “Gas Decompression Program” (GDP code).
4. Numerical analysis on rapid decompression process in base natural gas mixtures The proposed mathematical model is validated on the experimental data on rapid decompression in a shock tube having base natural gas mixtures (Botros et al., 2003). Rich gas mixtures are not considered in this study, because a liquid condensation from gaseous phase is usually occurred at rich and medium rich gases. The presented model does not account for a liquid condensation phenomenon in this study. The model will be extended into the case of two-phase gas–liquid flows in pipes, where a liquid condensation will be taken into account, and the model will be validated on the experimental measurements of rapid decompression in pipe systems having rich gas mixtures (Botros et al., 2007) in other work. Results of this work will be presented in future papers.
The experimental measurements of the decompression wave speed in natural gases were conducted by TCPL (TransCanada PipeLines) at TCPL Gas Dynamic Test Facility in Didsbury, Alberta, Canada (Botros et al., 2003). The main test section of the facility is a shock tube having a length of 30 m and the inner diameter of 49.325 mm. The internal surface of the shock tube has a roughness better than 40 micro-inches (i.e. less than 1.016 × 10−6 m). Basing on the value of the roughness-to-diameter ratio (ε/Dpipe = 20.5 × 10−6 ; ε = 1.016 × 10−6 m; Dpipe = 49.325 × 10−3 m), the internal surface of the shock tube is determined as the “smooth” tube surface (Botros et al., 2010b). Therefore, the wall roughness does not influence on the flow parameters strongly. The classical Blasius correlation (4)–(6) for the smooth pipes is allowed in this case. The influence of the internal surface roughness will be examined numerically in details by using different friction correlations, where the internal roughness will be taken into account, in other future work as well as other experimental measurements (Botros et al., 2010b) will be employed for the model validation. A rupture disc is placed at one end of the shock tube, which is upon rupturing. A decompression wave propagates up into the pressurized test section. A few high frequency responses pressure transducers (PT) are mounted into the tube in order to capture the time history values of the expansion fan (Botros et al., 2003). Several rupture discs were introduced into the shock tube end for different pressures at rupture. Decompression wave speed values are determined from the time between signals from P1 to P8 pressure transducers as well as the time between signals from P5 to P6 transducers. Fig. 1 shows the schematic view of the experimental decompression tube. Despite the fact that the main test section of the shock tube has a length of 30 m, the computational decompression pipe having a length of 50 m is simulated from the point of view of the numerical scheme stability. The inner diameter of the computational pipe is 49.325 mm. One end of the shock tube is selected to be the closed end. The rupture disc boundary condition is modelled at other pipe end. Distances between the rupture disc and pressure transducers are selected to be identical to real distances on the experimental
Table 1 – Initial pressure (MPa), temperature (K) and gas composition (mole %). Cases
Pinitial
Tinitial
N2
CO2
C1
C2
C3
i-C4
n-C4
i-C5
n-C5
Case 1 Case 2 Case 3
10.41 13.8 20.67
274.07 264.77 264.72
0.697 0.647 0.699
1.097 1.197 1.279
92.955 93.02 92.757
4.076 3.876 4.075
0.8 0.909 0.861
0.099 0.108 0.103
0.137 0.158 0.146
0.066 0.059 0.053
0.073 0.026 0.027
67
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
Pressure (MPa)
9 8 7 6 5
14
PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
12
Pressure (MPa)
PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
10
10
8
6
4 0
2
4
6
8 10
40
50
60
70
80
90
100
Time (ms)
decompression tube as it is shown in Fig. 1. Those distances between PT P1, P8, P5, P6 and the rupture disc are 0.4, 1.15, 16.4, 18.4 m, correspondently. Three experimental measurement data cases (Botros et al., 2003) having different initial values are selected in order to be simulated (Table 1). The following gas mixture composition (base natural gases), initial pressure and initial temperature are selected (Botros et al., 2003). Predictions are started with the initial pressure of 10.58 MPa and the temperature of 247.4 K in each computational cell of the pipe for the case 1 (Table 1). New values of the velocity, temperature, density and pressure are re-calculated after each time step. The upper limit of the time step is selected from the point of view of the numerical stability of calculations. The decompression is very quick process and basic parameters of gas mixture are changed very rapidly in time. Therefore, the time step has to be small enough in order to simulate continues change in pressure and temperature values in time at every computational cell. The computational time step is equal to 2 × 10−5 s in all of considered cases. Pressure values are collected at P1, P8, P5 and P6 PT locations after each time step. Fig. 2 shows the pressure time history values at P1, P8, P5, P6 PT locations for the case 1. Experimental points are shown in all figures as symbols. Continues lines represent predictions, which are made by using the proposed model. Fig. 3 shows measured (Botros et al., 2003)
1.0 Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
Pressure ratio (-)
0.8 0.7 0.6 0.5 50
100
150
200
250
5
10
15
40
60
80
100
120
140
Time (ms)
Fig. 2 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 1.
0.9
0
300
350
400
Decompression Wave Speed (m/s) Fig. 3 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 1.
Fig. 4 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 2.
and predicted values of the decompression wave speed, which are determined from P1 to P8 and P5 to P6 PT locations, correspondently. Pressure values are normalized on the initial pressure value before rupturing. Figs. 2–3 show a good agreement between experimental data and predictions, which are made by using GDP code. The proposed model simulates the slope in the pressure time history values correctly at each PT location and shows a good agreement with the measured data. The pressure shows rapid decrease at P1 and P8 PT locations. Those PT are located near the rupture end of the pipe. The friction does not influence on the decompression wave speed strongly here and, therefore, it does not slow it down significantly. The slope of the pressure time history curves at P5 and P6 PT locations is more flatter compared to P1 and P8 due to strong influence of the gas to wall friction on gas mixture release from the pipe. Calculations of the decompression wave speed in base natural gas mixtures in considered experimental conditions (Table 1; Cases 1–3) were made by TCPL (Botros et al., 2003) by using the analytical GASDECOM (Eiber et al., 1993) code. In order to compare those predictions with GDP calculations, GASDECOM’s calculated values are taken from Botros et al. (2003) and are put into Figs. 3, 5 and 7 of the paper. Broken curves in those figures represent GASDECOM’s numerical data (Botros et al., 2003). GASDECOM’s predictions are in good agreement with the proposed model and experimental decompression wave speed values, which are determined from P1 and P8 PT locations. However, GASDECOM does not able to predict the decompression wave speed values, which are determined from P5 and P6 PT locations, because the friction force is not taken into account in GASDECOM’s analysis. The proposed model (i.e. GDP code) shows ability in calculating the decompression wave speed, which is determined from P5 to P6 PT locations, correctly and with high order of agreement with the measurements. Other predictions of the decompression wave speed values are made by using the same shock tube topology (Fig. 1), but having much higher initial pressure before rupture started (Table 1, Case 2). The initial temperature and gas composition are up-dated accordingly to the experimental values, which were selected for the simulations (Table 1). Pressure time history values show a good agreement with experimental data in this case (Figs. 4–5), also. Predictions, which are made by using GDP code, show much better agreement
68
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
1.0
1.0 Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
0.8
Pressure ratio (-)
Pressure ratio (-)
0.9
0.9
0.7 0.6
0.7 0.6 0.5 0.4
0.5 0.4
0.8
Measured (P1&P8) GDP code (P1&P8) Measured (P5&P6) GDP code (P5&P6) GASDECOM
50
100
150
200
250
300
350
0.3 100
400
200
300
400
500
Decompression Wave Speed (m/s)
Decompression Wave Speed (m/s) Fig. 5 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 2.
Fig. 7 – Measured and predicted decompression wave speed values as a function of dimensionless pressure ratio for the case 3.
with the experimental data (Fig. 5) compare to GASDECOM’s calculations, which are taken from Botros et al. (2003) again. The last set of predictions of the decompression wave speed are made for the same shock tube topology (Fig. 1) but having the highest initial pressure value (Table 1, Case 3). The agreement between predicted pressure time history values and experimental measurements has a high order of magnitude (Fig. 6). The proposed model simulates the decompression process in base natural gas mixtures correctly and provides very detailed information about local distribution of basic flow parameters (spatial and time history, both), which is very helpful in engineering analysis. The prediction, which is made by using the proposed model, shows much better agreement with the experimental data compare to analytical GASDECOM model (Fig. 7). The analytical GASDECOM (Eiber et al., 1993) is one of the most frequently used engineering software in gas field applications performing quick simulations of the decompression wave speed in natural gases. However, this model does not account for the friction between the gas and the pipe’s wall. The flow behaviour in a pipe’s area far away from the rupture place is totally different compared to the rupture end area due to friction force contribution. It was observed experimentally, that the decompression wave speed, which is determined from P5 to P6 PT locations, farther away from the rupture end,
is much lower than the corresponding one (Figs. 2–7), which is determined from a closer pair of P1 to P8 PT locations (Botros et al., 2003). Therefore, the analytical GASDECOM program itself does not simulate the decompression process in natural gas mixtures correctly in the case of small-diameter shock tubes, where the friction between the gas and the wall plays an important role. GASDECOM’s values are over-predicted in this case. Its application area is limited on the large-diameter pipeline cases due to this fact. Simulations on rapid decompression process in dry and base gas mixtures were made by TCPL (Botros et al., 2007) by using the well-know and frequently used commercial 1D software OLGA, which is developed by SPT-group (Bendiksen et al., 1991). All predictions, which were made by using OLGA, are found in bad agreement with the experimental data. The frontal wave speed values, which are predicted by OLGA code, are significantly lower than the measured one (Botros et al., 2007). Thus, the proposed mathematical model of transient compressible thermal multi-component gas mixture flows in pipes predicts the decompression process in a shock tube having base natural gases much better than analytical and mathematical models, which are available from the open source literature. Moreover, simulations, which are made by using the presented mathematical model, are quick in time and allow simulating of large-scale industrial pipes and shock tubes.
22 PT-P1 PT-P1 PT-P8 PT-P8 PT-P5 PT-P5 PT-P6 PT-P6
20
Pressure (MPa)
18 16 14 12 10 8 6 0 2 4 6 8 10
30 40
50 60 70
80 90
Time (ms) Fig. 6 – Measured and predicted pressure time history values at P1, P8, P5 and P6 PT locations for the case 3.
5.
Conclusions
A one-dimensional transient mathematical model of compressible thermal multi-component gas mixture flows in pipes is presented in the paper. The set of mass, momentum and enthalpy conservation equations for gas phase is solved in the model. Thermo-physical properties of multi-component gas mixture are calculated by solving the Equation of State (EOS) model. The Soave–Redlich–Kwong (SRK–EOS) model is chosen. Gas mixture viscosity is calculated on the basis of the Lee–Gonzales–Eakin (LGE) correlation. Numerical analysis on rapid decompression process in base natural gas mixtures, which is performed by using the proposed mathematical model, is presented in the paper. The model is successfully validated on the experimental data on rapid decompression in a shock tube having base natural gas mixtures from. Three experimental cases with high, medium
chemical engineering research and design 9 1 ( 2 0 1 3 ) 63–69
and low initial pressure values before rupturing are simulated. The proposed mathematical model shows a very good agreement with the experiments in a wide range of pressure values. It has been shown, that the mathematical model predicts the decompression in base natural gases much better than analytical and mathematical models, which are available from the open source literature. The presented model is highly necessary and useful in pipeline design and flow assurance. The minimum of fracture arrest toughness of the pipe wall material may be determined on the basis of the Battelle two-curve method with taking into account of the proposed model together with fracture propagation speed model. The influence of the pressure, temperature, fluid composition, and pipeline diameter may be examined by using the presented model.
Acknowledgement The author would like to acknowledge the support of PETROSOFT-DC for giving possibility to use the computer program (GDP code) performing the mathematical modelling and numerical analysis on rapid decompression in a shock tube having base natural gases.
References Bendiksen, K.H., Maines, D., Moe, R., Nuland, S., 1991. The dynamic two-fluid model OLGA: theory and application. SPE Prod. Eng. V6 (N2), 171–180. Blasius, P.R.H., 1913. Das Aehnlichkeitsgesetz bei Reibungsvorgangen in Fluessigkeiten. 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. In: American Gas Association Proceedings (AGA-2003).
69
Botros, K.K., Studzinski, W., Geerligs, J., Glover, A., 2004. Determination of decompression wave speed in rich gas mixtures. Can. J. Chem. Eng. 82, 880–891. 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. Pres. Ves. Pip. 84, 358–367. 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. Pres. Ves. Technol. 132 (051303), 1–15. Botros, K.K., Geerligs, J., Rothwell, B., Carlson, L., Fletcher, L., Venton, P., 2010b. 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. Pres. Ves. Pip. 87, 681–695. Burlutskiy, E., 2012. Mathematical modelling of non-isothermal multi-component fluid flow in pipes applying to rapid gas decompression in rich and base natural gases. In: Proceedings of the International Conference on Fluid Mechanics, Heat Transfer and Thermodynamics–ICFMHTT-2012, vol. 61, Zurich, Switzerland, 15–17 January, pp. 156–167. Eiber, R.J., Bubenik, T.A., Maxey, W.A., 1993. GASDECOM, Computer Code for the Calculation of Gas Decompression Speed. Fracture Control for Natural Gas Pipelines. PRCI Report No. L51691. Eiber, R.J., Carlson, L., Leis, B., 2004. Fracture control requirements for gas transmission pipelines. In: Proceedings of the Fourth International Conference on Pipeline Technology, p. 437. Lee, A.L., Gonzalez, M.N., Eakin, B.E., 1966. The viscosity of natural gases. J. Petrol. Technol., 997–1000. Patankar, S., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York. Soave, G., 1979. Equilibrium constants from a modified Redlich–Kwong equation of state. Chem. Eng. Sci. 27, 1197–1203. Wallis, G.B., 1969. One-dimensional Two-phase Flows. McGraw Hill, New York.