富友 牛雪天 敏政

(兰州理工大学 能源与动力工程学院，甘肃 兰州 730050)

在输油、输水管道系统中，由于油泵、水泵机组及阀门众多，存在发生事故断电、阀门间启闭不协调的可能，这将导致管路系统形成复杂的瞬变流动现象。这种流动现象发生、发展的过程会导致整个系统瞬态过程中水力特性的改变，从而可能导致管道的爆破或塌陷、水力元件的严重破坏、输油/水管线大范围泄漏，甚至发生爆炸与火灾。鉴于这些严重的后果，需从机理上了解这一瞬变现象的本质即弥合性水击。当水击压力波呈现与流动相反方向传播时，会出现水头压力的降低，当下降至饱和蒸汽压附近时，介质会发生汽化，形成汽、液两相流。一旦形成的气(汽)相发生压缩溃灭，会产生明显的局部射流高压，这种伴随气(汽)相的发展与溃灭的水击现象被称为弥合性水击。

Baltzer[1]、Adamkowski等[2]分别对弥合性水击现象进行了理论分析及实验拍摄，分析中发现这种现象从物理过程上归属于空化范畴，但其发生、发展的过程更为复杂：由于液体始终处于被动拉伸状态，其空穴生长、压缩溃灭过程的时间更为短暂，同时在压力作用下空泡被压缩溃灭后形成射流所造成的破坏性更大。Tanaka等[3]发现：当弥合性水击出现时，压力上升阶段的次级压力波峰值将比首级峰值更高。Adamkowski等[4]发现，在弥合性水击与非弥合性水击的过渡区域，会产生一个时间极短的压力脉冲，这种现象出现时的边界条件较为特殊，时间较短，极易被忽视。随着研究的进一步深入，发现弥合性水击空穴压缩溃灭时所产生的压力激升，并非在第二个压力周期初期就呈现，而是在高压平稳末端形成更高的射流高压[5]，该压力峰值会超过Joukowsky理论计算值[6]，这也正是射流水击现象。对于这种射流流动现象的详细研究，是全面揭示弥合性水击内部机理及流动规律中急需解决的问题。

为此，本文采用双流体模型，通过对弛豫过程的相间作用机理进行分析，提出了一种改进的、可用于描述泡粒变化的弛豫过程计算方法，通过气液两相激波管算例验证该计算方法的稳定性及准确性，并针对弥合性水击的内部流动过程进行了分析，揭示了弥合性水击中的射流-缓冲机理。

1 数学计算模型

本文的数学计算模型采用气-液两相双流体模型进行构建。双流体模型由一个空穴的输运方程、两个连续性方程、两个动量方程、两个能量方程组成[7- 9]。模型中认为双相均为可压缩介质，并认为气相、液相分别具有不同的压力、密度。空穴输运方程(1)的存在保障了方程组呈现双曲型特性，模型方程如下：

(1)

(2)

(3)

(4)

式中：k=1代表气相，k=2代表液相；I代表液相与气相的过渡区域；α为空穴率；p为压力；ρ为密度；单位体积下的总能量E=e+u2/2，e为内能，u为速度；μ为压力弛豫系数。

1.1 相间参数的构造

由于两相之间的密度差异较大，因此需要一个相间区域(相间边界)来进行过渡。如图1所示，这是一个使气体和液体处于平衡状态的区域。该区域的相间参数可以用如下方程来表述[8,10- 11]：

(5)

(6)

1.2 状态方程

状态方程建立了压力、密度和能量之间的关系。该方程使用了加强气体方程的修正形式[10,12]：

(7)

(8)

Tref为环境温度，psat为饱和蒸汽压力。

图1 气液相间结构的示意图

2 数值计算方法

2.1 改进的Godunov- HLL数值计算方法

把方程(1)-(4)整理成如下形式，忽略空穴输运方程中由于弛豫现象产生的源项：

(9)

(10)

式中，U=(α1ρ1,α2ρ2,ρu,α1ρ1E1,α2ρ2E2)，F(U)=(α1ρ1u,α2ρ2u,ρu2+α1p1+α2p2,(α1ρ1E1+α1p1)u,(α2ρ2E2+α2p2)u)，Q(U)=(0,0,0,pIU,-pIU)。

在方程求解中，用到了HLL近似黎曼求解器。HLL近似黎曼求解器的数值表达式如下：

(11)

式中：FL、FR、UL、UR分别代表网格节点i(i=±1/2)左、右两端的边界值；S+和S-分别代表不同的波速取值，其表达式为[7,13]

(12)

数值计算中的差分格式的构建方法如下：

(13)

(14)

(15)

由于空穴输运方程不能写成守恒形式，因此本文采用改进的Godunov方案进行构建：

(16)

2.2 边界构建方法

(17)

对全微分方程进行积分并取一阶近似，可以得到如下表达式：

T-1(UP-UA)Δt

(18)

(19)

式中，P、A、B的取值如图2所示。

图2 特征值及特征线在空间坐标的分布趋势

进一步将方程(18)、(19)改写为如下形式：

T-1(UP-UA)Δt=T-1DAΔt

(20)

3条特征线上数值采用以下形式进行插值计算：

(21)

式中：UB代表当特征波速为负值时，向上游传递的信息源节点；UA代表当特征波速为正值时，向下游传递的信息源节点；UC代表A、B间网格上的数值节点。

2 弛豫计算模型

本文通过结合有限弛豫法和无限弛豫法的计算过程[14- 15]，推导出一种能够计算泡粒变化的双流体相间弛豫构建方法[16]。根据空泡静力学原理，当压力降低至饱和蒸汽压力附近时，流动区域并不一定发生空化，仅当液体压力降至psat-2S/R(S为表面张力，R为泡粒半径)以下时，空泡才可能开始继续增长。通常在室温下，直径为10 μm的空泡可以承受14 kPa的张力[17]。因此，结合泡动力学相关特性，将-2S/R引入弛豫构建方法中。由于弛豫过程与空间变量无关，它仅是时间的函数[7]，因此，可以将原始方程式

(22)

(23)

(24)

变换为如下形式：

(25)

(26)

(27)

式中，上标*表示弛豫变量，上标0表示计算变量。

根据空泡静力学原理[7,17]，气相和液相压力之差可以表述为

(28)

将方程进一步整理，可以得到：

(29)

(30)

(31)

3 验证算例

3.1 水气两相激波管

为了研究弛豫模型及数值算法的稳定性及准确性，选用水气两相激波管算例[7]进行验证。算例中，水气两相激波管被分成两部分，管内相对长度0～0.7范围内充满水，管内相对长度0.7～1.0范围内充满空气。由于双流体模型的数学特性，需要默认两侧分别存在体积分数在10-8以下量级的第二相。水气两相激波管的具体参数如表1所示。数值计算解与解析解在t=2.29 ms时的对比如图3所示，数值计算中CFL值定为0.8。从图中可以看出，在t=2.29 ms时，激波管内的混合密度、压力、速度分布与解析解具有良好的一致性。根据混合密度在0.8～0.9相对位置处台阶式分布匹配关系可以看出，该弛豫模型及数值算法具有良好的稳定性及准确性。

表1 水气两相激波管参数

3.2 弥合性水击实验分析-辛普森水击实验

为验证计算模型对带有相间质量变化物理过程计算的有效性，选用典型的弥合性水击实验-辛普森水击实验[5]进行对比分析。辛普森水击实验主要由4部分组成：上游水箱、管道、球阀和下游水箱。管道长度为36 m，内径为0.019 05 m，壁厚为0.001 588 m。管道在12.5 m位置有一个90°的弯头。

在下游水箱附近的36 m位置处存在一处球阀，球阀可以通过手动操作实现快速关闭。整个管道的坡度为1 m，沿管道有3个压力变送器，分别位于

图3 水气两相激波管数值解与解析解对比

9、27和36 m处，系统中波速为1 280 m/s。实验管道系统如图4所示，实验基本参数如下：流量0.332 m/s，上游水头34.54 m，水头损失0.33 m，摩擦系数0.031 5，理论水头变化43.332 m，关阀时间0.036 08 s。状态方程中的修正热力学参数见表2，Q为单位质量的能量。

通过表2的修正热力学参数来表达方程组中的状态方程。整个瞬变过程通过尾部球阀处的信息输入来计算整个管路中的瞬变特征，其方法为：首先在尾阀处测量压力变化，然后根据式(32)计算速度边界条件。当计算的速度变为0时，速度边界将保持为0[5]。

图4 实验系统示意图

表2 状态方程中的修正热力学参数[12]

(32)

从表2可知，实验中阀门关闭过程所消耗的时间小于压力波传播回阀门所需的时间(2L/c)，因此可以将其归类为直接水击实验。由于水击造成的压力变化值大于上游边界压力，因此它是具有水柱分离功能的直接水击。

弛豫模型与辛普森水击实验之间的水头变化比较如图5所示，从图中可知，计算值与实验值具有很好的一致性。低压(0.15 s)后的二次和三次压力波的响应时间和幅度都具有良好的精度。数值计算经网格无关性验证后，考虑计算成本，使用了720个网格单元，计算时间步长为10-7s。

3.3 弥合性水击射流-缓冲机理分析

由茹科夫斯基直接水击计算公式可知，当发生直接水击时，压力的升高值是一定的，即等同于第一个水击压力激升的数值60.32 m。从图5中可以看出，低压后的二次和三次压力波的数值远高于60.32 m。为阐明这一流动现象的内部机理，取不同时刻下的管路速度变化与压力变化特征进行对比分析。为了减少管路高程及管路摩擦损失对流动特性的影响，将辛普森水击实验中的高程、管路摩擦损失忽略，并将管路简化为无量纲的相对长度，0代表上游起点，1代表下游终止位置。不同时间节点下的管路速度与压力变化趋势如图6、7所示。

图5 辛普森水击实验水头变化信息

图6 管路的内流速度变化趋势

图7 管路的压力变化趋势

经分析发现，对于产生的二次压力高于理论计算值的弥合性水击而言，其流动规律可以分为直接水击阶段(0.01～0.08 s)、缓冲水击阶段(0.09～0.16 s)、射流水击阶段(0.17～0.24 s)。

1)在直接水击阶段下，速度的变化趋势为：先下降至0，再呈现反向的流动(速度为负值)，如图 6(a)所示。从图7(a)可知：在0.01～0.05 s期间，由于边界处的流体突然停止运动，管路流体速度逐渐降低，动能转化为压力能，压力则随时间的延长而升高，使液体被压缩，形成一股压缩波向上游进行传播；在0.06～0.08 s时，由于压差作用，流动开始自下游至上游反向流动，此时形成膨胀波向下游传递，使压力恢复。

2)在缓冲水击阶段的0.09～0.11 s时段，速度开始反向减小，逐渐接近0，但速度变化并未继续遵从直接水击阶段中速度恢复至0的流动特征，而是仍具有由下游向上游的方向(负值)，如图6(b)所示。这一流动变化的产生，是由于此时空穴泡粒在管路末端处不断析出，形成反向缓冲作用。从图7(b)可以看出，相对长度在0.4～1.0内，管路内流体绝对压力降低至饱和蒸汽压力以下，验证了空穴泡粒存在并起到缓冲反向作用的说法。

在0.12～0.14 s时段，流动趋势开始变化为由上游向下游进行，速度变化为正值。在相对长度为0.9～1.0的范围内，0.14 s时的速度出现较大的梯度变化，此区域伴随发生了空穴泡粒的压缩溃灭，但此时速度并未恢复至初始边界速度(0.332 m/s)，而同时期的压力变化信息显示，此区域压力由低压逐渐恢复至初始边界压力，验证了存在空穴泡粒压缩溃灭及其对流动起到缓冲作用的说法。

在0.15～0.16 s时段，速度相继出现较大的梯度变化，且出现两种速度变化叠加传播：由下游向上游的增速波和由下游到上游的减速波。其中增速波的传播要早于减速波。这种现象是由于管路末端的空穴泡粒发生压缩溃灭，造成整个管道中的流体产生瞬态加速，形成增速波由下游向上游逐渐传递。而当边界处的流体受到阻碍突然停止运动，压强升高使液体被压缩，形成升压波、减速波向上游传播。由图7(b)可以看出，同时期存在明显的压力升高，但其升高值约为60.32 m，远小于直接水击所产生的压力值77.87 m，因此将这一区间定义为水击缓冲区。

3)由图6(c)可以看出，在射流水击阶段的0.17～0.20 s时段，速度出现两种变化的叠加传播：由上游向下游的增速波和由上游到下游的减速波。其中增速波的传播要早于减速波。上一阶段的增速波在上游边界处反馈后形成反向传播(即上游向下游传播)。与此同时，压强在管内逐渐升高，由于上下游存在较大压差(边界压力为34.54 m)，导致流体反向流动，形成减速波由上游传播至下游。而由于增速波对边界的冲击形成了射流作用，导致压力的进一步上升，在0.20 s时形成了90 m左右的二次升高值(如图7(c)所示)，这一数值的出现，影响了采用直接水击方法计算压力升高值的准确度，因此这也是弥合性水击中最应该注意的部分。

在0.21～0.24 s时段，开始重复0.09～0.12 s 时段的内部流动变化规律，但由于二次射流压力波的存在，同时忽略了流动中的能量损耗，形成了与直接水击压力波不同的压力表现规律。

4 结论

本文采用双流体模型，通过对弛豫过程的相间作用机理进行分析，提出了一种改进的可用于描述泡粒体积变化的弛豫模型，采用气液两相激波管算例验证了该计算方法的准确性，并使用该方法对弥合性水击中的射流-缓冲现象进行了内部流动分析，得到以下结论：

(1)对于低压过后所产生的二次压力高于理论计算值的弥合性水击而言，其流动按规律可以分为直接水击阶段、缓冲水击阶段和射流水击阶段。

(2)对于缓冲水击阶段，由于空穴泡粒的生成，降低了初次水击对边界处的压力冲击，使得水击呈现被缓冲状态，压力升高值小于理论计算值，此时射流作用消失。

(3)对于射流水击阶段，空穴溃灭时产生的射流加速并未直接作用于初次水击对边界处的压力冲击，而是经过反向传递后，使得本应恢复至初始速度值的速度获得了一个加速作用，这一加速过程产生明显的波动向水击波传播的反方向传播，在初次水击基础上造成二次压力升高，此升高值要高于理论计算值。