一.   前言

瞬变流摩阻与许多要素有关，例如时间、频率、流量、压力以及管壁弹性等，其处理方法要比恒定流复杂得多。在传统的水击分析计算中，为了便于计算，通常采用稳态时的摩阻系(因)数来计算整个瞬变过程的摩阻，并且采用一阶近似的方法使摩阻和流量的关系线性化。对于摩阻起重要作用的系统，各泵站的压头大部分甚至全部都消耗在管道摩阻损失上(长距离的输油管道)，传统的近似方法会给计算结果带来较大的误差，而且由于误差的累加作用，可能根本得不到正确的解。

到目前为止，对瞬变流摩阻损失的机理和计算方法的研究还很不充分。虽然传统水击分析方法不能很好地模拟水击衰减过程，但却能够有效预测最大的水击升压和降压值。因此，寻找更为完善的瞬变流摩阻计算方法，进一步提高水击分析的精度，对于管道的安全运行具有重要的指导意义。

二.   传统瞬变流摩阻计算方法的改进

1   动态摩阻系数法

在我国石油储运工程界，常采用列宾宗公式进行摩阻计算^[1]：

式中h_f——水力摩阻损失；

f——列宾宗摩阻系数；

Q——管道流量；

m——流态指数；

l——管道长度。

传统水力瞬变分析通常采用稳定流状态下的f值来进行摩阻计算，误差较大。本研究根据瞬变状态下流量的变化情况，采用不同的公式来计算f，即f的分区计算。文献[1]给出了层流及紊流各区的计算公式。通常情况下，雷诺数(Re)是用来划分流态区的。对于一条具体的管道而言，管径与管壁的相对粗糙度均可视为定值，因此可将按Re划分流态区转化为按Q来划分。

层流区：

紊流光滑区：

紊流混合摩擦区：

紊流完全粗糙区：

式中ν——液体的运动粘度；

D——管道内径；

K——管道的当量粗糙度。

图 1给出了采用动态摩阻系数法进行水击分析的基本计算流程。

图  1  动态态摩阻系数法的计算流程

下载: 全尺寸图片 幻灯片

文献[1]指出，只要计算出的流态区域与实际相符合，用式(1)计算瞬态摩阻是准确的，或者基本上是准确的。

2   动态摩阻因数法

在紊流混摩区中，难以用式(1)的形式来计算摩阻损失，不得不求助于近似方法^[1]，这样势必会引入附加的误差。

计算紊流混摩区最准确的公式是国际上公认的科尔布鲁克(Colebrook)式^{[1, 2]}：

式中.：λ——达西摩阻因数。

文献[1]认为，紊流区的水力计算不分区，都采用式(6)，以利于工程计算，并可提高混摩区的计算精度。文献[2]指出，式(6)能够很好地符合Moody图，用其计算整个紊流区的λ值是切实可行的。将式(6)变形，建立Q与摩阻损失的关系，推导过程如下。

达西-韦斯巴赫(Darcy—Weisbach)公式表示为^[1]

由式(7)可以得到：

将式(8)变形为：

式中V——液流流速;

g——重力加速度；

s——单位管长的压头损失，即水力坡度；

a₁——系数。

将雷诺数纂和式(9)代入式(6)，于是有：

令，可将式(10)转化为：

式(11)建立起了流量与摩阻损失(水力坡度)的关系。已知摩阻损失可以直接计算流量，而已知流量计算摩阻损失则需要迭代。图 2给出了采用摩阻迭代法进行水击计算的基本流程。虽然式(6)的解算过程稍显复杂，但可忽略分区问题，计算精度也能得到保证。

图  2  瞬变摩阻迭代计算流程

下载: 全尺寸图片 幻灯片

三.   辛普生方法处理摩阻积分项

水击基本微分方程中的摩阻项通常有两种表达形式，即和，现采用后者进行分析。根据辛普生方法^{[3, 4]}，令：

式中时刻的流速值。

将式(12)代入正向和负向特征线差分方程，得到：

式中A——管道截面积;

a——水击波速;

Q_L——特征线AP中点的流量；

Q_R——特征线BP中点的流量(见图 3)。

图  3  差分分计算网格

下载: 全尺寸图片 幻灯片

对Q_L和Q_R的求解不采用辛普生方法，因为特征线AL、LC及RC、RB中点的流量始终是未知量，这样是无法求出Q_L和Q_R的。考虑到计算方便和保证精度的要求，这里采用将摩阻项处理成线性隐式的方法^{[1, 5]}。

对于特征线AL：

对于特征线CL：

对Q_R的处理方法类似。采用这种方法，计算精度比一阶近似模型要高，而且算法也比较容易实现。求出Q_L和Q_R后，再联立式(13)和式(14)，可解出H_P和Q_P。

辛普生方法的截断误差为：

其中t₁≤ξ≤t₂辛普生方法计算精度高的原因除了其本身具有三阶精度外，还引入了Q_L和Q_R两个参数。这表明，在计算点P的物理量时，不仅考虑了前一时步相邻两点A、B的影响，还考虑了前一时步同一点流动状态的影响，这与水击波传播的实际情况是相符的。而采用一阶近似模型时，完全没有考虑AP或BP间流量的变化，精度必然很差。采用二阶精度模型时，虽然考虑了流量的变化，但是同一点的物理量，要经过两个时步之后才能产生影响，这与真实的物理过程有偏差，其计算精度自然也会受到一定影响。

在与传统方法相同计算时步下，采用辛普生方法，水击计算精度虽明显优于传统方法，但其计算量会相应增加，这主要是由于需要计算中点流量的缘故。如果再加上瞬态摩阻的计算时间，则总的计算量还将迸一步增加。考虑到本文方法的计算精度比传统方法高，因此，即使在较大的计算时步下，此方法仍可获得与传统方法相同甚至更高的精度。对于长输管道水击计算而言，采用该方法可以适当加大步长，在保证计算精度的同时，增加的计算量可能并不多甚至没有增加。因为计算量与管道分段数N的平方成正比，增大步长可以显著地减小计算量^[6]。

四.   算例

图 4为一简单管道，管道内径D=0.15 m，上游为水库，下游接球阀，稳态时球阀全开，自由流出，波速a=1 000 m/s，计算终端阀门关闭时，管道末端的最大水击压头。

图  4  简单单管道不惹图

下载: 全尺寸图片 幻灯片

表 1列出了管道长3 000 m时三种方法的计算结果，表中的n表示管道的差分段数。根据不同初始流量和上游压头，分别用三种方法计算了管道末端的最大水击压力，每种方法又分别采用不同的计算步长进行了计算。

表  1  管道(L=3 000 m)终端阀门处的最大水击压头

下载: 导出CSV 

| 显示表格

从表 1中的数据可以看出，变摩阻系数法或动态摩阻因数法取n=3时，其计算结果接近于传统方法取n=20时的结果；而取n=10时，又十分接近于传统方法取n=50时的计算结果。这说明为了得到相同的计算精度，用变摩阻系数法或动态摩阻因数法，再结合辛普生方法，可以将计算时步加大为传统方法的5~7倍左右，从而显著节省计算时间。

表 2列出了不同长度的长输管道终端关阀时，末端的最大水击压头。从表 2数据可以看出，在保证与传统方法具有相同计算精度的前提下，应用动态摩阻因数法计算所需要的n值远小于传统方法，并且管道越长，所需n值越小。这说明，对长距离管道而言，应用该方法能够大大减小计算工作量，同时还能保证相同甚至更高的计算精度。

表  2  管道末端的最大水击压头

下载: 导出CSV 

| 显示表格

五.   结论

从两个方面对传统瞬变流摩阻计算进行了改进，即在计算水击基本微分方程中的摩阻项时，紧密结合流量变化，采用动态摩阻系(因)数，并结合辛普生法计算瞬变摩阻，提高了计算精度，进而提高了长输管道水击数值分析的准确度。计算单个节点的压力和流量时，采用该方法的计算工作量要大于传统方法。但对于整条管道，为了获得相同计算精度的解，传统方法所需的计算节点数是该方法的数倍，对长距离管道而言，甚至高达数十倍，算例中的计算结果也充分说明了这一点。因此，该方法可以用较大的时间步长来进行管道水击分析，其计算工作量不但没有增加，甚至还会大大降低，这在长距离管道水击计算上体现得更为明显。