CC BY
10
УДК 532.517+532.542
Вестник СПбГУ. Сер. 10. 2016. Вып. 4
N. N. Ermolaeva
COMPUTER MODELLING OF THE SEA GAS-PIPELINE GLACIATION AND OF THE FLOW CHARACTERISTICS BEHAVIOR IN UNSTEADY REGIMES
St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation
The mathematical model of non-stationary gas mixture transmission by the sea gas pipeline is suggested taking into consideration the possibility of pipeline glaciation. The algorithm of the numerical solution of a nonlinear model equation system is presented. The algorithm has been implemented in C+—+ in the form of program complex, which enables to calculate time changes of the temperature, density, pressure and speed of gas mixture in every section of the gas pipeline and time changes of an ice layer thickness on outer surface of the gas pipeline. As an example, the calculation of one of variant gas transmission is illustrated by the chart, which shows variations of the ice layer thickness for different sections. Refs 6. Figs 5.
Keywords: sea gas-pipelines, gas transmission, dynamic of glaciation, nonstationary flow, mathematical models, Lax—Wendroff-type scheme.
Н. Н. Ермолаева
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ОЛЕДЕНЕНИЯ МОРСКОГО ГАЗОПРОВОДА И ПОВЕДЕНИЯ ХАРАКТЕРИСТИК ПОТОКА В НЕУСТАНОВИВШИХСЯ РЕЖИМАХ
Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Предложена математическая модель нестационарной транспортировки смеси газов по морским газопроводам, учитывающая возможное оледенение газопровода. Приведен алгоритм численного решения нелинейной системы уравнений модели. Он реализован на языке С+—+ в виде комплекса программ SGPIT, позволяющего рассчитать изменения во времени в каждом сечении газопровода температуры, плотности, давления и скорости потока газовой смеси, а также толщину слоя льда на внешней поверхности газопровода. Пример расчета одного из вариантов транспортировки газа проиллюстрирован графиками, демонстрирующими изменения толщины слоя льда во времени для разных сечений. Библиогр. 6 назв. Ил. 5.
Ключевые слова: морские газопроводы, транспортировка газа, динамика оледенения, нестационарное течение, математические модели, схема Лакса—Вендроффа.
Introduction. The Arctic seas of Russia contain big natural gas fields. Work has continued on the development of gas condensate fields in these seas. One of the problems is transportation of the extracted gas from an offshore platform to the mainland. A computer model of flow of a real multicomponent gas mixture in a sea hyper-pressure gaspipeline and at possible glaciation of the outer surface of pipeline is an effective instrument to analyze and to compute the sea gas pipeline parameters and the gas transmission regimes. Computer simulation of these processes plays an important role in the long-term
Ermolaeva Nadeczda Nikolaevna — PhD of physical and mathematical sciences, associate professor; n.ermolaeva@spbu.ru
Ермолаева Надежда Николаевна — кандидат физико-математических наук, доцент; n.ermolaeva@spbu.ru
© Санкт-Петербургский государственный университет, 2016
forecasting of gas pipeline work in the steady and the unsteady regimes. The mathematical models of the steady-state regimes for these problems are presented in the book [1]. The present work continues these studies for unsteady regimes of the gas transportation through the sea gas pipeline operating under such conditions that glaciation of outer surface of the pipeline is possible.
Mathematical model of gas flow and of gas pipeline glaciation dynamics (model 1). The model 1 described of the formulas
+ + ^ = (3,
dt ' dz\r \ pj) R
e = e + v2/2, (4)
-- kpT CP'2 (5)
1 ~5p (1 +Sp)VT' 3 c
e = cvT--j-j=Hl + 5p), (6)
the initial conditions for p, T, u, (7)
the boundary conditions for p, T, u; (8)
the unit of calculation q, when no ice layer (unit A)
dT
-^- = a1C{T1), r G (i?, i?i), (9)
r = R: T\ = T, q = Ai^i,
or
r> rp rp dT2
r = R\ . Ii = l2, Ai—— =A2^—, or or
dT2 dt
n rp rp ^ dT2 , dTw
r — K 2- 12 — J-w, — —,
or or
dTw
=aw£(Tw), r e (R.2, i?2 +
dt
r = R2 + S* : Tw = T*, the initial conditions for Ti, T2, Tw ; the unit of calculation q, when ice layer exists (unit B) the equations and the conditions (9)-(12) are as in unit A,
10) 11) 12)
13)
14)
15)
16)
r = R2: T2 = Th А2^ = АЗ, (17)
or or
3T-
гб(Д2,Д2+е), (18)
r = R + С : Ti = Tw = T*, (19)
dTw
= Л/ Па ■
dt
А ^
i dr
— —
^2+« dr
1йЛ, (20)
Й2+«
dT
-^ = aw£(Tw), re(R2+^R2+^ + S„), (21)
dt
r = R2 + e + : Tw = T*, (22)
the initial conditions for Ti, T2, Ti, Tw. (23)
In presented model equations (1)-(3) are equations of the continuity, of the momentum and of the energy, (5) is the Redlich—Kwong equation of state, (6) is the caloric equation [2]; (9), (11), (14), (18), (21) are the heat equations in the pipeline wall layers, in the thermal boundary layer and in the ice layer, respectively; (20) is the Stefan condition.
In system (1)-(23) we use the following designations: u,p,p,T are the flow velocity, the density, the pressure, and the temperature of a gas mixture, correspondingly, which are functions of time t and coordinate z coinciding with the gas pipeline axis; e = e(z,t), e = e(z, t) are the mass densities of energy and internal energy; A is the hydraulic resistance coefficient; q is the radial component of flux density of the internal energy (the heat flux vector) upon the inner surface of pipeline in z-th cross-sections; h, c, S are the dimensional constants in the Redlich—Kwong equation of state determined by a given chemical composition of a gas mixture [3]; cv in equation (6) is the specific heat of an ideal gas (including ideal gas mixtures); r is the radial coordinate in cylindrical coordinate system (r,(fi,z); C = ^^ (r^) is the Laplace operator in cylindrical coordinate system (r, z, ip) for ^ = ^ = 0; £(t) is the ice thickness in z-th cross-sections; au = Ak/(pk Ck), Ak, pk, ck u Tk = Tk(r,t) are the thermal diffusivity, the thermal conductivity density, the specific heat and the temperature distribution in k-th layer, indexes k = 1, 2, w, i correspond to the following regions: 1 to the first layer of pipeline wall (steel), 2 to the second layer of pipeline wall (concrete), i to the ice layer, w to the thermal boundary layer of water; 7 is the latent heat (of fusion), R is the inner radius of the gas-pipeline, R1, R2 are the outer radiuses of the pipeline wall layers; T* is the ambient temperature, T* is the seawater-ice transition temperature; S*, S** are the thicknesses of thermal boundary-layer when no ice layer and when one exists, respectively.
The thermal boundary-layer thickness depends on the many factors, particularly, on bottom currents and on ground coupling. And generally, the thickness of the thermal boundary-layer changes when the outer wall pipeline glaciation occurs. It is difficult enough to calculate theoretically the values of thicknesses S*, S**. This information can be obtained from the external flow problem of pipeline.
The estimation of these values in steady-state regimes is presented in book [1]. In practice following inequalities hold: S** < S* C R2. Model 1 can be extended to a larger number of layers.
In model 1 the processes in gas flow are described in terms the averaged (over the cross-sectional area) pressure, in the velocity, in the temperature and etc. For model 1
it is assumed that the heat exchange processes with the surroundings and the ice formation processes have parametric dependence on coordinate 2 (through dependence the gas temperature T(z,t) on 2). This presumption is justified, if the radial component of the temperature gradient in pipeline wall layers much greater than the axial component.
Conditions of the glaciation beginning. During solving general non-steady problem for every z can be calculated the time moment t(z) in which in this section the ice layer occurs. For this two conditions must be satisfied
Л2 —
dr
T2R2 Л) < (24)
dTw
> К
R2,i dr
■ (25)
R2,i
Inequality (25) follows from the Stefan condition (20) at £ ^ 0 for such pipeline cross-sections that was no the ice layer in initial time. For the pipeline cross-sections in which is no the ice layer the nonsteady regimes are calculated using system (1)-(8) and unit A. In other cross-sections the system (1)-(8) is added by unit B.
Solution algorithm. Numerical solution algorithm for system (1)-(23) is based on the below Lax—Wendroff scheme.
At each time step the heat flux q is calculated either by the equations system of unit A or of unit B. Equations of unit A are solved by an implicit method using elimination method. To numerically equations solution of unit B is used a method with the explicit tracking of moving surface [4], which is an iterative finite difference method, with variable time steps. In this approach the time step size is variable and, it is determined so that the ice thickness increased on given constant value during this time step.
As demonstrated previously [5] and as confirmed by calculations of test cases using model 1, in studied problems, when inlet gas temperature is significantly higher than the ambient temperature, the length of pipeline l* exists, for which in segment z < l* the intense heat exchange with the ambient occurs.
In this segment the substitution of unit A by the quasi-steady state approximation can cause errors in the calculations. For z > l* the calculation of unit A can be made using quasi-steady heat transfer description with sufficient precision.
For an illustration in Fig. 1 have plotted the temperature distribution T(z) along the pipeline rote at time t = 10 hour. The temperature is calculated using model 1 for the unsteady treatment of unit A (the curve 1) and for the quasi-steady model (the curve 2). From Fig. 1 is obvious that beginning l* ^ 250 km the unit A calculation can replace by its quasi-steady version.
The pipeline glaciation occurs at the end of the pipeline route, in the region z >> l* where is possible to use the quasi-steady approximation to calculation the temperature in the pipeline wall layers. Besides, the glaciation rate of outer pipeline surface, as expected, is much less than the temperature change rate in wall layers. Both above factors allow to replace the unit B by its quasi-steady version.
The quasi-steady version of unit B:
Tk(r)= Ak + Bk lnr, k = 1, 2,i,w, (26)
T(z,t) = Ai + Bi ln R, (27)
A1B1 = A2B2 = XiBi, (28)
A2 + B2 ln R2 = Ai + Bi ln R2, (29)
T, K
Fig. 1. Gas temperature distribution T(z) along the pipeline rote at time t =10 hour: for the unsteady heat transfer model (1 ) and the quasi-steady model (2 )
Ai + Bi lnR + £) = Aw + Bw lnR + £) = T*, À^I^ — ÀWBW _
r~2Tl
T* = Aw + Bw ln(R2 + £ + S**) The coefficients A& (t), Bk (t) satisfying to system (26)-(32) are equal:
T* - T(z,t)
Bi =
In — + — In — + — In ^ ^
R À2 Ri Xi R2
, Bi = (Xi/Xi )Bi,
Ai = T(z,t) - Bi lnR, B2 = (Ài/À2)Bi, Ai = T* - Bi ln(R2 + £),
A2 = T* + Bi
Bw =
(T* - T* )
R2 + £ + S**
ln
R2 + £
The expression for Bw (34) can be simplified by use
S*
Bw —
R2 + £ (T* - T* )(R2 + £)
< 1:
(30)
(31)
(32)
(33)
(34)
**
Stefan condition (31) turns into following ODE for £(t):
, Bj (T*-n)
^Тё"*™-= (35)
The desired expression for heat flux q in pipeline sections, in which ice layer exists, in the quasi-steady approximation is given by the equation
AiBi
1=—Д-- (36)
Here B1 is calculated from equation (33) and ice thickness C(t), which is a member of B1, is calculated from ODE (35).
As noted above, in most practical cases the inequalities 6** < 6* с R2 hold. For calculation of the glaciation beginning time t(z) is permissible to use the quasi-steady versions of units A and of B for sections z > Z*. In this case, the conditions (24), (25) can be represented in the form of one inequality
T(z, t) < T* - (T* - T.) (f + ^ fin § + ± In f ) ) . (37)
\6** Ai 6** \ R Л2 Ri J J
The time t(z) for every cross-section is determined as a moment in which the gas temperature T(x,t) satisfies the condition (37) in this cross-section. (In considered problems, the gas temperature is monotonically decreasing function in every cross-section.)
Boundary and initial conditions. We consider non-stationary problem of gas trans-portation, in which nonstationarity is due to the gas consumption variations and due to the ice formation processes.
Initial conditions. For model 1 the initial conditions are the distributions of the density p0(z), of the temperature T0(z) distributions and of the ice thickness £o(z) in steady state regime:
Q
t = 0 : у = pu = const = ——г,
nR2
p(z) = P0(z), T (z) = T0(z), £,(z)= Co (z).
The functions p0(z), T0(z), £0(z) are calculated using steady variant of model 1 presented in the book [1], Q is the mass flow rate, which is constant for the steady regime.
Boundary conditions. The gas flow in the gas pipeline is subsonic. In considered problem the unchanged over time the inlet pressure and the inlet gas temperature are given. Using these values the density and the internal energy are determined from the caloric equation and the equation of state. At outlet the law of variation of the specific flow rate y*(t) is given. Thus, the boundary conditions are written this way:
z = 0: p(0,t) = p0, e(0,t) = £0,
z = L : y(L,t) = y* (t),
L is the gas pipeline length.
The algorithm of the numerical solution of model 1. Turn from the velocity u to the flow rate y = pu and cast equations (1)-(6) in non-dimensional form in the variables p, y, p, £ using for the dimensionless quantities the same notation:
dt dz '
ду д (у2 \ у2
m + d~z{~ + mip) = "m2 7'
д ( y2\ д ( y3 y
— [ре + тз — j + — I ye + m3 -2 + mA-p ] = m5 q,
p = Ш6
dz PT
Ш7
P
(1 - miop) e(T, p) = mgT - mg
(1 + m10p)VT ' ln(l + mwp)
VT ■
The dimensionless complexes m1 — m10 are expressed through the physical parameters and the characteristic values from the formulas
m1
m-4 =
Px
pxuXx
m-2
Px£x
m5 =
4R 2Ai txTx
px^xRrx
Ш3
2£x
m6 =
hpxTx
m-7
mg
2
CPx
PxVTx
3c
mg
mio = 5px,
25exy/T.
lx, tx are the characteristic length and time, px, px, Tx are the characteristic density, the pressure and the temperature of the real gas mixture, ex, ux are the characteristic internal
energy and the velocity of gas mixture (e
xx
3 c
cvTx--j-j=r\n{l + 5px), ux = Q/(pxTvR2));
the quantities px,px,Tx are connected by the equation of state (5).
For the numerical solution of model 1 we use modified Lax-Wendroff scheme [6], which appeared to be the most preferable for the considered problems by the count rate and the simplicity of implementation.
The algorithm consists of two steps. At every step the desired values of the density, of the flow rate and of the internal energy are determined explicitly:
stage I:
stage II:
where
Un+1/2 _ 0 5 (Un + Un ) г?'
Uk+l/2 U'° \Uk + uk+1) | tk+1 -
0.5t
Д
k+1/2,
Tjn+1 rrn Fn+1/^ Fn+1/2 Uk ~ Uk . fe+1/2 Pk-1/2
Г Д
n+1/2
U=
p y
y
pe + m3 — \ p
F =
y
y\
--h mip
p 3
у , У
ye + + mA-p
p2 p
0
2
-m2
\ m5 q /
2
u
x
p
x
p
x
Cv T x
£
x
Here n, t are the number and the time step size, A, k are the number and the space grid step size 2. The temperature and the pressure are found from the preceding equations (5), (6).
In this scheme the calculation ql+i/2 is carried out according to equation (10) or (36). As follows from formula (36), the calculation x'n+i/2 includes the ice thickness xl+i/2 at the (n + 1/2)-th time step. This value is determined using equation (35), which can be written in form
d£ si (T* - T(z,t))
dt (R2+0 D +
R2
- S2 = f(T(z,t),C), (38)
here £, RX are the dimensionless quantities, but the temperature T,T*,T* was left dimensional
Ai tx txXw (T * — T)
«1 = -5-, «2 = -7-,
rX YPi rxYPiO**
D — — In — + — In —
Ai R AX Ri
The desired ice layer tickness is found explisit from the discrete counterparts of
equation (38).
The Lax—Wendroff is a second-order difference method in both time and space.
This algorithm of model 1 solution has been implemented in C++ in the form of programm complex SGPIT. It allows to calculate all flow characteristics in any time and ice layer thickness on the outer surface of the gas pipeline. It allows to calculate all flow characteristics in any time moment and the ice layer thickness on the outer surface of the gas pipeline by given parameters of the pipeline construction.
As an illustration we selected the boundary condition variant, for which the inlet pressure is equal to 17.2 MPa to demonstrate the glaciation dynamics. It should be pointed out that by the inlet pressure more then 21 MPa, while the other parameters remained the same, the segment of glaciation and the ice layer thicknesses are small.
The calculation result using program complex SCRIPT for the calculation unsteady non-isothermal flow of gas mixture under glaciation conditions is presented at the following values of the parameters:
Q = 570 kg/sec; T* = 272.15 K; T* =271.15 K; L = 450 km; R = 0.5 m, Ri = 0.54 m, RX = 0.66 m, S* = 0.020 m, S** = 0.018 m;
Ai = 24 W/(m-K), Ax = 1.7 W/(m-K), Aw =0.56 W/(m-K), Ai = 2.3 W/(m-K);
Y = 335000 J/kg, ci = 450 J/kg-K, cX = 924 J/kg-K, cw =4200 J/kg-K; pi = 928kg/m3, pi = 10000kg/m3, pX = 2300kg/m3, pw = 1005kg/m3;
T(0,t) = 303.15 K, p(0,t) = 17.2 MPa, A = 0.0083, c.v = 1712.25 J/kg-K.
In calculations we used following dimensional values of the characteristic parameters: tx = 3600 sec, rx = 0.01 m, px = 138.02 kg/m3, Tx = 283.15 K, lx = 10 km, them correspond the characteristic pressure px = 15.2 MPa and the characteristic velocity ux = 3.69 m/sec. The values h = 502.9, c = 12 297.58, S = 0.0018 were chosen correspond to the gas mixture with more methane. The calculations were carried out by the dimensionless time step t = 0.000035 and by the dimensionless space step a = 0.01. In Fig. 2 the dependence time
t from the cross-section coordinate 2 is shown. t is the time of beginning glaciation in the part of gas pipeline, in which at the initial time t0 the ice layer is absent. By the abscissa axis, 2 in km is indicated, by the ordinate axis time in hours is indicated. In Fig. 3 the temperature change of gas T along the pipeline route at t = 5 days is shown.
hour
T, К
Fig. 2. Dependence time t from the cross-section coordinate z
400 500
z, km
Fig. 3. Temperature change of gas T along the pipeline route at t = 5 days
In Fig. 4 the change of ice layer thickness £ along the pipeline route at t = 5 days is shown. In Fig. 5 the change of ice layer thickness in cross-sections l = 350 km (see a) and l = 450 km (see b) during five days.
С cm
0 100 200 300
Fig. 4. Change of ice layer thickness £ along the pipeline route at t = 5 days
400 500
z, km
с cm
а
b
0.7
4
0.6
0.5
0.4
0.3
0.2
0.1
0
О 20 40 60 80 100 120 140
О 20 40 60 80 100 120 140
t, hour
Fig. 5. Change of ice layer thickness £ in the cross-sections l = 350 km (a) and l = 450 km (b) during first five days
By the abscissa axis, time in hours is indicated, by the ordinate axis ice layer thickness in centimeters is indicated.
Conclusion. The mathematical model of the unsteady non-isothermal turbulent gas flow on the sea gas pipeline including the conjugate model of glaciation dynamics of outer surface pipeline was presented.
The effective algorithm of numerical solution of model 1, implemented in the form of complex program SGPIT, was suggested.
Results analysis of computer simulation of the gas transportation processes on the sea gas-pipeline and the conjugate processes of its outer surface glaciation was introduced.
The region of flow, in which can use the quasi-steady approximation of the heat transfer model was determined.
For the considered problems the admissibility of quasi-steady approximation by the simulation of glaciation dynamics of outer gas-pipeline surface was proved.
References
1. Kurbatova G. I., Popova E. A., Filippov B. V. Modeli morskih gazoprovodov [Models of sea gas-pipelines]. Saint Petersburg, Saint Petersburg University Publ., 2005, 156 p. (In Russian)
2. Ermolaeva N. N., Kurbatova G. I. Analiz podhodov k modelirovaniju termodinamicheskih processovv v gazah pri visokih davlenijah [The analysis of the approaches to the modeling of thermodynamic processes in gas flow at hyperpressure]. Vestnik of Saint Petersburg University. Series 10. Applied Mathematics. Computer science. Control processes, 2013, issue 1, pp. 35-45. (In Russian)
3. Reid C., Prausnitz M., Sherwood K. The properties of gases and liquids. New York, Saint Louis, San Francisco, MeGraw—Hill Book Company Publ., 1977, 688 p.
4. Vasilev F. P. O metode konechnih raznostej dlja reshenija odnofaznoj zadachi Stefana [On finite difference methods for the solution of Stefan's single-phase problem]. Zh. Vychisl. Mat. Mat. Fiz., 1963, vol. 3, no. 5, pp. 861-873. (In Russian)
5. Ermolaeva N. N., Kurbatova G. I. Kvasiodnomernaja nestacionarnaja model processov v morskih gasoprovodah [Quasione-dimensional non-stationary model of processes in a sea gas-pipeline]. Vestnik of Saint Petersburg University. Series 10. Applied Mathematics. Computer science. Control processes, 2015, issue 3, pp. 55-66. (In Russian)
6. Roach P. J. Computational Fluid Dynamics. Albuquerque, Hermosa Publ., 1976, 618 p. (Russ ed.: Roach P. J. V'ichislitelnaja gidrodinamika. Moskow, Mir Publ., 1980, 616 p.)
For citation: Ermolaeva N. N. Computer modelling of the sea GAS-pipeline glaciation and of the flow characteristics behavior in unsteady regimes. Vestnik of Saint Petersburg University. Series 10. Applied mathematics. Computer science. Control processes, 2016, issue 4, pp. 75—85. DOI: 10.21638/11701/spbu10.2016.407
Статья рекомендована к печати проф. Н. В. Егоровым. Статья поступила в редакцию 5 июня 2016 г. Статья принята к печати 29 сентября 2016 г.
