基于动态摩阻模型的大落差输油管道不稳定流动模拟

李旺，冯亮, 黄立宇, 廖代兵, 陈莎, 贾文龙, 孙璐璐

(1.国家管网集团西南管道有限责任公司，成都 610000；2.西南石油大学石油与天然气工程学院，成都 610000)

中国西南地区的中缅原油管道、兰成渝成品油管道、兰成原油管道都是大落差输油管道，管道高低点的局部落差高达500～1 500 m，压力差值高达4～8 MPa[1]。特别是在管道停泵、关阀等不稳定流动工况下，管道局部高点处容易出现低压，使油品中的轻质组分挥发，从而在高点处出现气泡聚集，形成液柱分离。过低的压力容易导致管道压瘪，而气泡破裂后容易形成弥合水击导致管道内出现超压和爆管。准确分析大落差输油管道在不稳定工况下的压力、流量变化是合理预判弥合水击现象的基础和前提。

在管道不稳定流动模拟[2]中，水力摩阻系数是计算管道压力及流量损耗的重要参数。传统的摩阻处理方法多采用拟稳态摩阻假设，即管壁切应力采用Darcy-Weisbach摩阻公式计算[3]。然而，在管道的不稳定流动过程中，油品在管道内始终处于不断加速、减速状态，流体的惯性会和管壁之间产生相比于稳态流动更大的摩阻损失(即动态摩阻)[4]。实验和理论研究表明，忽略不稳定流动中流体惯性影响(即动态摩阻)，会低估水击摩阻值，不能准确预测水击压力波的衰减过程和波形畸变[5-7]。

动态摩阻系数与多因素有关，计算较为复杂，许多学者对此进行了实验和理论研究[8]，并提出了多种模型。Zielke[9]基于Navier-Stokes方程，提出了层流瞬变流的管壁切应力计算模型，该模型计算简单，不依赖经验参数，对计算与频率相关的摩阻损失精度较高，但需要考虑每一时步的流速值，计算时间长、内存消耗大。Trikha[10]采用加权函数改进了Zielke模型，使得计算只需前一步的值，简化了计算过程。Martin等[11]将管道分为近壁层流区和核心区，推导出了适用于水力光滑管中任意雷诺数湍流非恒定摩阻计算模型。Brunone等[12]、Vardy等[13-14]考虑瞬时当地加速度和对流加速度，得到了非恒定摩阻计算模型，更真实地模拟了不稳定流摩阻。近年来，有学者注意到了动态摩阻对水利工程的影响[15-17]，研究了动态摩阻对输水管道泄漏过程压力和流量畸变的影响[18-20]，但国内外对于动态摩阻对输油管道的不稳定流动过程的影响研究相对较少。

基于以上研究，在管道不稳定流动连续性方程、动量方程和能量方程基础上，现引入Brunone-Vitkovsky动态摩阻模型[21]，结合特征线方法[22]和Newton-Raphson方法，对大落差输油管道不稳定流动工况进行模拟，分析大落差管道水力瞬变特征。

1 数学模型

原油在长输管道内的流动可视为一维流动。根据流体流动的质量守恒、动量守恒、和能量守恒原理，建立管道内的一维不稳定流动模型为

(1)

式(1)中：U为自变量列向量，U=[P,V,T]T；A为系数矩阵；B为列向量。稳态流动条件下，原油单位面积管内壁的摩擦力τ可按Darcy-Weisbach[13]摩阻公式计算。对于不稳定流动，τ按Brunone-Vitkovsky动态摩阻公式计算[9-11]，公式为

(2)

A=

(3)

B=

(4)

式中：kt、kx分别为当地加速度和迁移加速度系数；sign(V)为符号函数，当V≥0时，sign(V)=1；当V<0时，sign(V)=-1；ρ为介质密度，kg/m3；P为管内压力，Pa；θ为管道与水平面间的倾角，rad；t为时间变量，s；V为介质流速，m/s；d为管道内径，m；g为重力加速度，m/s2；a为流体波速，m/s；τ为流体与单位面积管内壁的摩擦力，Pa；x为管道位置变量，m；K为管道的总热传递系数，W/(m2·K)；Cv为介质比热容，J/(kg·K)；T为介质温度，K；T0为地温，K；f为水力摩阻系数；D为管道外径，m。

2 模型求解

2.1 特征线方程

基本方程[式(1)]为拟线性双曲型偏微分方程组，它有3条特征线，可以沿特征线将其转化为全微分方程，然后对全微分方程积分便可得易于数值处理的全微分方程。从物理意义上讲，特征线即为压力波和温度波的传播路径。为此，建立了基于特征线法的模型求解格式。根据线性代数理论，系数矩阵A的3个特征值表达式为

(5)

式(5)中：

(6)

(7)

如考虑kx=kt=k3，且a>>V，则有

(8)

(9)

k3按式(10)取值为

(10)

式(10)中：对于层流，C=0.004 76；对于紊流，有

(11)

由Li(A-λiE)=0可得每个λi的相应特征向量Li为

(12)

由此可得3个特征值对应特征线方程及其相容性方程如式(13)～式(15)所示。

前向特征方程C+：

(13)

后向特征方程C-：

(14)

将连续性方程带入温度特征方程，可得温度特征方程C为

(15)

2.2 有限差分方程

管道瞬变流动计算要求给出管道沿线的压力、流量和温度随时间变化的规律，为此将管道沿时间和长度方向离散成矩形网格，如图1所示。为保证解的稳定性，网格长度必须满足Courant条件[13]，即

图1 管道不稳定流动差分网格示意图

(16)

对于上述特征方程中的相容性方程沿各自的特征线进行积分，对摩阻项采用二阶近似来防止解的不稳定。可以得到C+、C-、C相对应的差分方程为

(17)

(18)

(19)

式中：Q为管道内输量，m3/s；A为管道截面积，m2；下角标参数N、S、R、H见图1；其他参数见第1节。

的特征线与上一时间层网格线的交点。S、N和R点的值要通过对上一时步的值进行线性插值而得到。上述方程仅能用于计算管道内节点的参数，并且以稳态时的水力参数作为初始值。对于上游边界点，需要将C-上的特征线方程与边界条件联解得到。对于下游边界，需要将C+和C上的特征线方程与下游边界条件联解得到。对于中间泵站、阀门等非管元件，需要将C+、C-、C上的特征线方程与泵的特性方程，阀门方程联立求解得到。上述由边界条件、以及有限差分方程构成的非线性方程组可采用Newton-Raphson迭代法进行求解。

3 算例分析

3.1 算例描述

某大落差加热输油管道全长438 km，管道沿线高程及水热力分布图如图2所示。管道规格为Φ273.1 mm×8.7 mm，不设保温层，管道总传热系数取2.6 W/(m2·K)，环境温度为30 ℃，全线共设置4个热泵站，各热泵站高程、里程位置如表1所示。管道最高点为P(100.94，3.424)km，沿线最大落差为549 m。稳定流动时管道内流量为418.6 m3/h；为减少油品密度、黏度等参数对于水力计算结果的影响，取输送油品相对密度为0.82，动力黏度为120 mPa·s。

表1 站场布置及进出站压头、温度参数

图2 管道沿线高程及水热力分布图

3.2 计算过程

基于考虑动态摩阻的输油管道不稳定流动模型，采用C#语言编制了输油管道不稳定流动计算程序，分析大落差管道的不稳定流动特征。设定管道起点、终点边界条件分别为压力温度和流量边界条件，结合泵站的进出站压头，计算得到沿线压力、温度和流量分布。以此稳态条件作为初始值，分析No.3泵站停运的瞬态工况。停泵运行工况具体设置如下：分析No.3泵站在1 s时两台泵同时停运，同时当No.3站上游压力大于下游时旁通阀门打开工况下的权限压力和流量变化。分别按照拟稳态摩阻和动态摩阻方法计算管道内的压力和流量变化，对比高点P与No.4站的进站压力及流量的差异，结果如图3～图6所示。

图3 高点P处的压力变化

图4 高点P处的流量变化

图5 No.4站进站压力变化

图6 No.4站进站流量变化

3.3 结果分析

图3和图4的计算结果表明，在No.3站两台泵同时停运后，其增压波经过87 s传递到最高点P(图3)，在此之前由两种方法计算P点的压力、流量一致，表明在稳定流动条件下动态摩阻对管道水力参数的影响可忽略不计；87 s后，由两种方法计算的P点压力、流量开始出现差异，最大差异出现在208 s时刻，此时P点波动最剧烈(流量、压力变化最快)，压力差0.024 4 MPa(图3)，流量差6.01 m3/h(图4)。基于稳态摩阻和动态摩阻计算的压力、流量之间的最大相对偏差分别为2.06%、2.11%，且基于稳态摩阻计算的压力、流量值偏高。

图5和图6的计算结果表明，No.3站在1 s后两台泵同时停运，其减压波经过162 s传递到No.4站，之前由两种方法计算No.4站进站压力、流量也是一致的。但是，162 s后，由两种方法计算的No.4站压力、流量开始出现差异，最大差异出现在300 s时刻，此时No.4站处波动也是最剧烈的(流量、压力变化最快)，压力差0.027 6 MPa(图5)，流量差1.49 m3/h(图6)。基于稳态摩阻和动态摩阻计算的压力、流量之间的最大相对偏差分别为3.24%、0.43%，且基于稳态摩阻计算的压力、流量值偏高。

模拟结果表明，当管道运行过程中发生不稳定流动时，由两种方法计算的压力和流量会出现较大的差异。这是因为不稳定流动过程中的动态摩阻是由拟稳态摩阻项和附加摩阻项之和进行表征的，动态摩阻计算方法考虑了附加摩阻项对管道不稳定流动过程的影响，使计算结果更为准确的反应管道不稳定流动过程，这与式(4)是吻合的。由于基于动态摩阻计算的管道压力低于稳态摩阻值，因此能够更加准确地预测在大落差管道高点可能出现的负压，进而分析其可能引发的弥合水击等现象，为采取更加切实可行的管道安全保障措施提供理论和技术支撑。

4 结论

(1)考虑输油管道不稳定流动过程中由流体不稳定流动带来的流体惯性影响，基于一维流动的质量、动量、能量守恒方程，以及由稳定流摩阻和附加项摩阻构成的Brunone-Vitkovsky动态摩阻模型，建立了大落差管道不稳定流动模型。

(2)采用特征线法、有限差分法，求解输油管道不稳定流动数学模型，分析了某大落差输油管道中间泵站停运后站前局部高点与下游泵站进站的压力、流量变化趋势，对比稳态摩阻模型与瞬态模型的结果差异。

(3)瞬态摩阻模型计算局部高点处的压力、流量较稳态摩阻模型的计算值之间的最大偏差分别为2.06%、2.11%，下游泵站进站压力、流量的最大偏差分别为3.24%、0.43%，且基于稳态摩阻计算的压力、流量值偏高。基于动态摩阻模型对大落差管道不稳定流动工况进行模拟，对于提高管道压力、流量等工艺参数的预测精度，制订更加有效的不稳定流动安全防护措施具有重要意义。