Arch. Min. Sci., Vol. 58 (2013), No 1, p. 131144 Electronic version (in color) of this paper is available: http://mining.archives.pl DOI 10.2478/amsc-2013-0009 ANDRZEJ OSIADACZ*1, FERDINAND E. UILHOORN*, MACIEJ CHACZYKOWSKI* ASSESSING HYDRATE FORMATION IN NATURAL GAS PIPELINES UNDER TRANSIENT OPERATION OCENA ZJAWISKA TWORZENIA SI HYDRATÓW W WARUNKACH NIEUSTALONEGO PRZEPLYWU GAZU W GAZOCIGACH This work presents a transient, non-isothermal compressible gas flow model that is combined with a hydrate phase equilibrium model. It enables, to determine whether hydrates could form under existing operating conditions in natural gas pipelines. In particular, to determine the time and location at which the natural gas enters the hydrate formation region. The gas flow is described by a set of partial differential equations resulting from the conservation of mass, momentum, and energy. Real gas effects are determined by the predictive Soave-Redlich-Kwong group contribution method. By means of statistical mechanics, the hydrate model is formulated combined with classical thermodynamics of phase equilibria for systems that contain water and both hydrate forming and non-hydrate forming gases as function of pressure, temperature, and gas composition. To demonstrate the applicability a case study is conducted. Keywords: hydrates, pipeline, natural gas, transient gas flow W artykule omówiono model nieustalonego, nieizotermicznego przeplywu gazu w rurocigu, który uwzgldnia model gazowego hydratu w stanie równowagi fazowej. To pozwala okreli czy hydraty mog tworzy si w okrelonych warunkach eksploatacji gazocigu a w szczególnoci okreli czas oraz miejsce ich tworzenia. Przeplyw gazu jest opisany za pomoc ukladu równa róniczkowych czstkowych utworzonych w oparciu o równanie zachowania masy, pdu, energii oraz równanie stanu wykorzystujce równanie Soave-Redlich-Kwonga. Za pomoc mechaniki statystycznej, model hydratu jest formulowany w oparciu o równowag fazow dla ukladów zawierajcych wod oraz gazy tworzce i nie tworzce hydraty jako funkcj cinienia, temperatury oraz skladu gazu. Slowa kluczowe: hydraty, rurocigi, gaz ziemny, nieustalony przeplyw gazu GAS ENGINEERING DIVISION, FACULTY OF ENVIRONMENTAL ENGINEERING, WARSAW UNIVERSITY OF TECHNOLOGY 00-653 WARSAW, NOWOWIEJSKA 20, POLAND CORRESPONDING AUTHOR. E-mail: email@example.com 1. Introduction Hydrate formation in natural gas pipelines should be avoided since they can plug the pipeline. Hence, hydrate control is fundamental in flow assurance. This article focuses on hydrate tracking in high-pressure pipelines under transient conditions. A non-isothermal, transient gas flow model and hydrate phase equilibrium model are solved. The gas pressure and temperature as function of time are calculated from the gas model and used to evaluate whether they are outside hydrate formation region. The hydrate model is based on statistical mechanics approach introduced by Van der Waals (1956) and Van der Waals and Platteeuw (1959). By deploying, the classical thermodynamics as boundary of hydrate kinetics, the pipeline operator is able to demarcate the region labeled as hydrate risk. Thus, irrespective of the kinetics, i.e., nucleation and growth, hydrates do not form, if the pressure and temperature of the gas are outside the hydrate stability region. The boundary is defined by the equilibrium hydrate formation curve. The formation point on the curve is the moment when nucleation starts and hydrates will start to form. The gas flow is described by a set of partial differential equations resulting from the conservation of mass, momentum, and energy. Although, many hydrate models have been published, none of the models is combined with transient gas flow modeling. Nomenclature Roman symbols a Equation of state mixture parameter, radius spherical core (Å), A Cross-sectional area, m2, as Isentropic wave speed, m/s, b Equation of state mixture parameter, ci[1, 2,3] Mathias-Copeman coefficients, Ck,m Langmuir constant of molecule k in cavity m (Pa-1), cp Specific heat at constant pressure, J/kg ·K, 0 cp Ideal isobaric heat capacity, J/kg ·K, d Diameter pipe, m, f Friction factor, -, fugacity, Pa, g Gravitational acceleration, m · s2, E g0 Excess Gibbs energy from UNIFAC model, J, H Specific enthalpy (J/mol), hydrate phase, I Power value of pressure and temperature, -, k Boltzmann constant, J/K, K Pipe roughness, m, L Liquid phase, m Pressure effect, -, · m Mass flow rate, kg/s, Rm Free radius of cavity type m (Å), Rs Specific heat constant, J/kg ·K, pc Critical pressure, Pa, q Heat flow into the pipe, J/m · s, qn Gas flow at normal condition, m3/s, q1 R Re Rm p ph pn qn t T Tc Ts U v V w yk, m z Mixing rule constant for PSRK, Gas constant, J/mol·K, Reynolds number, Free radius of cavity type m (Å), Pressure, Pa, Hydrate formation pressure, Pa, Pressure at normal conditions, Pa, Gas flow at normal conditions, m3/s, Time, s, Temperature, K, Critical temperature, K, Surrounding temperature, K, Overall heat transfer coefficient, W/m2 · K, Velocity, m/s, Molar volume, m3 · mol, Frictional force per unit length and time, N/m, Fractional occupancy of cavity m by guest molecule k, Compressibility factor, -, coordination number, -, Greek symbols Ice phase, Empty hydrate lattice, i Activity coefficient of component i, Depth of the potential well (J) JT Joule-Thomson coefficient, K/Pa,, Inclination angle of pipe, radian, Isentropic exponent, µ Dynamic viscosity, Pa · s, chemical potential (J/mol), Density, kg/m3, Cores distance at zero potential (Å), Number of cavities per water molecule, Spherically symmetrical potential energy (J). 2. Non-isothermal transient flow model From the laws of conservation of mass, momentum and energy the basic equations in term of partial differential equations describing a one-dimensional transient flow are expressed as follows (Thorley, 1987) ¶p ¶p a 2 ¶v +v + r as = s ¶t ¶x ¶ x cp T æ T æ ¶z æ æ æ q + wv æ ç1 + ç ç ç ç ç ç z è ¶ T èp ç è A è è è (1) ¶v ¶v 1 ¶p w +v + =- g sin (q ) ¶t ¶x r ¶ x Ar (2) ¶T ¶T a s2 æ T +v + ç1 + ¶t ¶x cp ç z è 2 æ ¶ z æ æ ¶ v as = ç çç è ¶ T è pç ¶x cp p è æ p æ ¶z æ æ æ q + wv æ ç1 - ç ç ç ç ç ç z è ¶ p èT ç è A è è è (3) where the isentropic wave speed is defined as as2 = zRsT é æ æ ê1 - p ç ¶z ç - p z è ¶p èT rcp T ê ë æ T ç1 ç z è æ ¶z æ æ ç ç ç è ¶ T è pç è 2ù (4) ú ú û and the frictional force per unit length is given by w= fr v v pd 8 (5) In this work, the friction factor f is calculated from Techo et al. (1965) as follows æ 1.964 ln (Re) - 3.8215 K æ f = -0.8685 ln ç + ç Re 3.71d è è -2 (6) As a matter of convenience Eqs. (1-3) can be written in terms of mass flow. This is accomplished by using the state equation for a real gas p = zRT r (7) The resulting set of equations is 2 ¶p as = ¶ t cp T . . . æ T æ ¶ z æ æ æ q mzRT æ é mzRT as2 m æ p æ ¶z æ æ ù ¶p wç - ê ç1 + ç ç ç ç + ç1 - ç ç ç ú 2 ç z è ¶T è p ç è A ç pA pA è z è ¶ p èT ç û ¶ x è ë pA ê è è èú . 2. æ 2 a m T æ ¶z æ æ ¶T a s ¶m - s ç1 + ç ç ç ç TA è z è ¶T èp ç ¶x A ¶x è (8) ¶T a2 = s ¶t cp p . . æ p æ ¶z æ æ æ q mzRT æ mzRT ¶T as2 æ T æ ¶ z æ æ ç ç1 - ç ç ç ç + w ç1 + çç ç pA2 è pA ¶x cp ç z ç ¶T çp ç è èè è z è ¶p èT è è A è . . . é mzR æ T æ ¶ z æ æ ¶T mTR z æ p æ ¶ z æ æ ¶p zTR ¶ m ù ´ê - 2 ç1 - ç ç ç + ú ç1 + ç ç ç ç z è ¶T èp ç ¶ x p A ç z è ¶p èT ç ¶ x pA ¶x ú ê pA è è è è ë û (9) . . m ¶m =¶t T . . æ T æ ¶ z æ æ ¶T m æ p æ ¶ z æ æ ¶p m 2zR æ T æ ¶ z æ æ ¶T + ç1 - ç ç ç ç1 + ç ç ç ç1 + ç ç ç ç z è ¶T èp ç ¶x z è ¶T èpç ¶t pç z è ¶p è T ç ¶t pA ç è è è è è è .2 . . æ m TRz æ pAg sin (q ) p æ ¶z æ æ æ ¶p mzTR ¶m + ç 2 ç1 - ç ç ç - Aç -wç ç ¶x ç pA ç z è ¶p è T è è pA ¶ x zTR è è (10) The heat transfer term q represents the amount of heat exchanged between the gas and the surroundings per unit length and per time and defined as follows q = dU(T Ts) (11) 3. Real gas effects The compressibility factor is calculated from the predictive Soave-Redlich-Kwong (PSRK) group contribution method (Holderbaum & Gmehling, 1991). This method uses the Soave-RedlichKwong (SRK) equation of state incorporated with the modified Huron-Vidal first-order mixing rule (Dahl & Michelsen, 1990). The UNIFAC model (Hansen et al., 1991) is used as the excess Gibbs energy for the mixing rule. The isobaric heat capacity, cp, and isentropic exponent, , are obtained numerically from the group contribution method as follows æ ¶p æ Tç ç V æ ¶2p æ è ¶T è V 0 cp = c p + T ç 2 ç dV -R æ ¶p æ ¥ è ¶T èV ç ç è ¶V è T (12) and é R æ æ ¶z æ æ ù k = ê1 çz +T ç ç ç ú è ¶T è p ç ú ê cp ç èû è ë -1 (13) The dynamic viscosity of the gas is a function of pressure and temperature and obtained from Katz et al. (1959). µ(p, t) = m · µ(pn, T ) (14) where m accounts for the pressure effect and µ(p0, T ) is the dynamic viscosity at 101.325 kPa. 4. Hydrate model By means of statistical mechanics, Van der Waals and Platteeuw (1959) derived the chemical potential of water in the hydrate phase, Dm w ( b -H ) æ æ (b ) (H ) = m w - m w = - RT ånm ln ç1 - å yk, mç è è m k (15) The phase equilibrium criteria, which are derived from entropy maximization, states that the chemical potential between the hydrate phase H and the coexisting water phase (liquid L or ice ), must be equal, Dm w (b - H ) ( b b H L = m w - m w = m w - m w or a = Dm w b -L or a) (16) where µw is the chemical potential of the metastable empty hydrate lattice. Instead of using (H (L the chemical potential as equilibrium criteria, the fugacity is deployed, i.e., f w ) = f w or ). The fugacity of water in the hydrate phase is defined as (H ) fw æ ææ æ (b ) = f w exp ç ånm ln ç1 - å yk, mçç èè è k è m (17) ( where f w ) is the fugacity of the empty hydrate lattice and yk, m the fractional occupancy. The latter is the product of Langmuir constant and fugacity, expressed as yk, m = 1 + å Ci, m fi Ck , m fk (18) The Langmuir constant accounts for the gas-water interactions in the cavities and is given by 4p Ck , m º kT Rm 0 æ w i, m ( r ) æ 2 exp çç r dr kT è è (19) It is directly related to the configurational partition function, which depends on the intermolecular potential and integral over the interaction volume. The spherical cell potential i, m(r) of component i in cavity m is obtained from é s 12 æ æ s6 æ a a 5 æù i wi, m ( r ) = 2 zme i ê 11 çd i10 + i d i11 ç - 5i çd i4m + i d i , mç ú ,m ,m , Rm Rm êRi r è è Rmr è èú ë û and (20) N i, m = 1 N -N -N éæ æ ai æ ai æ ù r r êç1 ç - ç1 + ç ú êè Rm Rm è è Rm Rm è ú ë û (21) where Rm is the cell radius of cavity m and z the co-ordination number of the cavity. N is 4, 5, 10 or 11. The Kihara parameters, i, i and ai of component i were fitted from experimental data summarized in Sloan (1998). The chemical potential of liquid water is formulated as p L æ ém (b,0) - m (w ,0)ù T éH (b ) - H (L)ù æ éV (b) - V (L)ù ( fw(L) = f wb ) exp ç ê w ú- ê ú dT + ê ú dp - ln ( xwgw )ç 2 çê ç RT0 RT û ú T0 ë RT p0 ë û û èë è (22) and for pure ice p a æ é m (b ,0) - m (w ,0)ù T éH (b ) - H (a)ù éV (b ) - V (a)ù æ (a (b fw ) = f w ) exp ç ê w ú- ê ú dT + ê ú dpç 2 çê RT0 RT û ç ú T0 ë RT p0 ë û û è èë (23) The activity coefficient is calculated from the UNIFAC model. The composition of the guest molecule in the liquid phase is obtained from Henry's law, xi, w = f i, g /Hi, w where f i, g is the fugacity of the gas component i and xi is the mole fraction of the gas dissolved in water. The effect of pressure on Henry's constant Hi, w is given by Krichevsky and Kasarnovsky (1935). (H) (L Now equating, f w = f w or ) and by iteration, the hydrate formation pressure or temperature for a given gas composition is calculated. A multiphase flash calculation is performed to calculate the phase of the mixture of known total composition. The PSRK group contribution method calculates the fugacity of all components in vapor and liquid phases. The activity coefficient i E and excess Gibbs energy g0 are calculated from the UNIFAC model. The procedure above assumes free water is present. In case of two-phase equilibria the algorithm is slightly modified whereas the hydrate pres(H) L L sure or temperature is calculated by equating, f w = f w,g where f w,g is the fugacity of water in ( the hydrocarbon calculated with the PSRK method. The fugacity of the empty hydrate, f w ) is obtained from multivariate regression. The correlation yields ln ( fw (b ) ) = å ai pI i =0 an i ,0 I i ,1 (24) where I is the power of pressure and temperature and ai the corresponding coefficient. The values are given in Table 1. TABLE 1 Coefficients of Eq. (24) Coeff. × 106 1 sI 2 sII 3 Ii, 0 4 Ii, 1 5 a1 a2 a3 a4 a5 112.2364 13.7154 10443.5035 2725049.9757 59305.0641 111.3305 12.5241 9765.9689 2619579.9389 59967.2289 a6 a7 a8 a9 a10 452.5014 237589591.4157 7962184.1570 113972.4068 726.6272 313.3178 235909838.5132 7749197.5651 166339.7300 3841.3551 The method of solution that needs to solve the set of non-linear equations uses a modification of Powell's (1970) hybrid algorithm. The n-point Gaussian-Legendre formula is used to avoid singularities and to achieve an efficient integration of the cell potential function. The model is coded in FORTRAN and verified with experimental data. The AADP for single and mixed gases in the presence of free water is 8.86%. For two-phase equilibria, the AADPPMV is 6.07%. The model has been simultaneously compared with CSMHYD, a thermodynamic model developed by Sloan (1998). Details of the hydrate model used in this paper, including verification with experimental data can be found in Osiadacz et al. (2009) and in Osiadacz et al. (2012). 5. Case study In the case study, the gas flow and hydrate model are combined to assess whether the gas pressure and temperature operate within the hydrate risk region. It enables to calculate the time and location of intersection with the hydrate formation curve. Several parameters are taken from the offshore Balgzand-Bacton pipeline. In the calculations, the following data are used: Gas: The gas is a mixture with a molar composition of CH4 = 98.3455, C2H6 = 0.6104, C3H8 = 0.1572, i-C4H10 = 0.0299, n-C4H10 = 0.0253, i-C5H12 = 0.0055, n-C5H12 = 0.0040, N2 = 0.0303 and CO2= 0.7918. The density n = 0.695 kg/m3 and the dynamic viscosity µ = 12.59 µPa · s. Pipe: The distance between the compressor stations is L = 235 km and the pipe diameter do = 1.016 mm. The roughness off the internal coating is assumed 15 µm. The simulations are carried out for three arbitrary selected heat transfer coefficients, which are presented in Table 2. TABLE 2 Overall heat transfer coefficients (Guo et al., 2005) Insulation type W/m2 · K Solid polypropylene Pipe-in-Pipe syntactic polyurethane foam Glass Syntactic Polyurethane The boundary conditions are p(0, t ) = 13.5 MPa T(0, t ) = 318.15°C q n (L, t ) = f (t ) (25) where f (t) is depicted in Figure 2 with a time interval, t [0,744 h]. The period function is taken from historical flows measured at the Balgzand-Bacton compressor station, i.e., the month March 2010. The initial values are obtained by setting the rates of change with time in Eqs. (1-3) equal to zero, i.e., v/t = 0, p/t = 0 and T/t = 0. The resulting equations for a horizontal pipe are ¶p = ¶x as2 cp T æ T ç1 + ç z è 2 æ ¶z æ æ æ q + wv æ a s w ç ¶T ç ç ç A ç + Av è èp ç è è è æ as2 æ çv - ç vè è (26) ¶T a2 = s ¶ x vcp p æ p æ ¶ z æ æ æ q + wv æ a s2 ç1 - ç ç ç ç çç z è ¶p èT ç è A è vcp è è æ T ç1 + ç z è 1 ¶p æ æ ¶z æ æ æ w ç ç ç çç è ¶T è p ç è vAr vr ¶ x è è (27) The ordinary differential equations are solved by the fourth-order Runge-Kutta method. The partial differential equations (8-10) are solved by the method of lines (Schiesser, 1991). A five-point, fourth-order finite difference approximation is used for the spatial derivatives and the stiff non-linear ODE system was solved numerically with a second-order accurate implicit Runge-Kutta algorithm. Fig. 1. Boundary condition at x = L A typical water specification lies between 64 and 112 mg/Nm3 (Carrol, 2003). The corresponding hydrate loci including free water are depicted in Figure 2. Fig. 2. Hydrate locus. (--) Free water, (-) 112 mg/Nm3 and (-- --) 64 mg/Nm3 -- To reduce the computation time, i.e., to find the intersection with the pressure and temperature of the gas and corresponding time and location, the results were fitted to the following equation: c æ æ ph (T ) = exp çc1 ln (T ) + 2 + c3 T + c4ç T è è where the constants are given in Table 3. (28) TABLE 3 Coefficients for Eq. (28) Free water 64 mg/Nm3 112 mg/Nm3 c1 c2 c3 c4 9382.6983 1340924.26 16.2874 53082.8267 3752.6514 535138.0822 6.6204 21219.082 5687.8204 808171.2124 10.0643 32133.8729 It should be noticed, that the gas flow model does not take into account the presence of free water. It is assumed that the small amount of water has a neglectable effect on gas pressure and temperature. Figure 3, depicts the moment and location at which the pressure and temperature in the gas pipeline enter the region of hydrate risk. It is the time and location when the pressure and temperature intersects with the hydrate formation curve (Fig. 2). Assuming free water in the pipeline and U = 2.84 W/m2 · K, the location at which the gas enters the hydrate region is ranging between 60 and 90 km along the gas pipeline. By improving the pipeline insulation down to U = 0.96 W/ m2 · K, the moment of entering the hydrate risk region starts at 173 km. The intersection curve is truncated at 100, 300 and 600 h. It means that the pressure and temperature are outside the hydrate risk region. If U is further reduced to 0.17 W/m2 · K, the gas remains outside the risk region. There is no intersection with the hydrate formation curve. Similar results are shown in Figure 4. However, reducing the water content to 112 mg/Nm3, does not significantly improve the situation. The results are similar as for free water in the pipeline (Figure 3). Figure 5, depicts the location as function of time for a pipeline with U = 2.84 W/m2 · K. Due to the lower water content of 64 mg/Nm3 and overall heat transfer coefficient, the natural gas remains outside the hydrate stability region. Fig. 3. L(t) at which the gas enters the hydrate region, assuming free water present in the natural gas. -- (--) U = 2.84 W/m2 · K and (-) U = 0.96 W/m2 · K Fig. 4. L(t) at which the gas enters the hydrate region, assuming a water content of 112 mg/Nm3. (--) U = 2.84 W/m2 · K and (-) U = 0.96 W/m2 · K -- The moment of intersection for free water, 112 mg/Nm3 and 64 mg/Nm3 with U = 2.84 W/m2 ·K is shown in Figure 6. As mentioned before, it can be clearly seen that the difference between free water and 112 mg/Nm3 is small. The solid polypropylene insulation with U = 2.84 W/m2 · K (Table 2) is for all specified water contents not sufficient to exclude the risk of hydrates. Figure 7, shows the results assuming U = 0.96 W/m2 · K and only in combination with a water content of 64 mg/Nm3 hydrates do not form a risk. Fig. 5. L(t) at which the gas enters the hydrate region, assuming a water content of 64 mg/Nm3 and U = 2.84 W/m2 · K Fig. 6. L(t) at which the gas enters the hydrate region with different water contents, assuming U = 2.84 W/m2. (--) Free water, (-) 112 mg/Nm3 and (-- --) 64 mg/Nm3 -- Fig. 7. L(t) at which the gas enters the hydrate region with different water contents, assuming U = 0.96 W/m2. (--) Free water and (-) 112 mg/m3 -- 6. Conclusion In this work, transient gas flow and hydrate modeling are combined in order to assess the risk of hydrate formation. A case study has been conducted for an offshore gas pipeline is order to demonstrate its application. Different water contents together with heat transfer coefficients have been assumed to find the time and location at which the pressure and temperature enter the hydrate risk region. By combining hydrate with gas flow modeling it is possible to determine time and location at which the natural gas enters the hydrate stability region. Moreover, to evaluate if the specified water content and pipeline insulation are sufficient to avoid the risk of hydrate formation. Further research is required to verify the assumption that the flow can be assumed as a single phase. In the gas flow model, the small amount of water in gas pipeline has been neglected when calculating the gas temperature and pressure.
Archives of Mining Sciences – de Gruyter
Published: Mar 1, 2013