液柱分离状态下水锤二次升压现象计算与研究

富 友 解朝蓬 蒋 劲

(1.兰州理工大学 能源与动力工程学院,兰州 730050;2.安徽三联泵业股份有限公司,安徽 马鞍山 243000;3.武汉大学 动力与机械学院,武汉 430072)

在液柱分离水锤现象中,二次升压现象一直是研究的重点内容.这种压力激升往往伴随着空泡的生成与溃灭,在阀门处尤为明显,其压力激升甚至比根据Joukowsky公式计算的最大压力还要高得多[1].有学者分析了二次升压过程中阀门处空泡的初生—发育—溃灭的过程,从而探究空泡溃灭对流动特性的影响[2]；也有试验结果表明,二次升压现象在低瞬变强度下会尤为突出[3],但都对其内部压力波的变化规律研究较少.

在伴有液柱分离的瞬变流动分析中,空泡的生成使得水锤压力波的变化程度与单相水锤不同.探索更加精确的瞬变流动摩阻模型对水锤压力波进行捕捉,有助于分析二次升压压力波的内部变化规律[4].此外,这种二次升压的产生会对机组、阀门和仪表产生较大的危害,所以对其进行研究是非常有必要的.本文基于双流体数学模型,采用数值模拟与试验对比相结合的研究方法,构建出一种改良的水锤数值解法,目的是为了探讨二次升压的内部流动机理,分析边界水头、边界压力、水锤波速、管道长度和管道直径对二次升压的影响.

1 计算方法

1.1 双流体数学模型

本文采用欧拉-欧拉双流体模型[5-7],该模型中两相均为可压缩.模型包含一个空泡输运方程、两个连续性方程、两个动量方程和两个能量方程.表达式如下：

方程中：下标1代表气相；下标2代表液相；α为空穴率；u为速度(m/s)；ρ为密度(kg/m3)；p为压力(Pa)；J u为非恒定摩阻；J s=fu|u|/2D；D为管道直径(m)；J g为重力项,J g=gsinθ；μ为压力弛豫系数；λ为速度弛豫系数；E=e+u2/2；pI为两相相间压力,pI=α1p1+α2p2；uI为两相相间速度,uI=(α1ρ1u1+α2ρ2u2)/(α1ρ1+α2ρ2)；α1+α2=1.

为使上述方程组封闭,引入状态方程：

其中：γk为比热容比,γk=C p,k/C v,k；C p,k为定压比热容；C v,k为定容比热容.π1通常为0,π2表达式见式(9)：

其中：Tref为试验环境温度；psat为饱和蒸汽压.

波速方程为：

方程(8～10)中：k=1代表气相；k=2代表液相.

1.1.1 相间弛豫系数

在使用双流体模型时,必须考虑瞬态流动中的压力和速度弛豫问题.其决定了模型的计算稳定性[8-9].相间弛豫系数包括速度弛豫系数和压力弛豫系数.

速度弛豫系数可表达为：

在压力弛豫系数中引入空泡静态平衡条件：

其中：S为表面张力；R为泡粒半径.

进而改进的压力弛豫系数表示为：

1.1.2 非恒定摩阻模型

瞬态流动模型可采用单系数和双系数摩阻模型来计算压力波的衰减.单系数和双系数摩阻模型由于时间和空间偏微分项的存在,使其在双流体模型中的应用变得复杂.加权函数非恒定摩阻模型[10]是对Zielke模型和Trikha模型[11]的简化.基于该模型的数值编码,使用加权函数,可以节省大量计算所需的时间.

对Trikha模型简化,通过加权平均,将瞬态流动阻力模型转化为如下形式：

式中：

其中：υ为运动粘度(m2/s)；g为重力加速度(m/s2)；y i为速度变量(m/s)；y i中i取1、2、3,对应m i和n i中3个数.

1.2 求解双流体模型的数值方法

1.2.1 Godunov-Rusanov格式

将方程组(1～7)进行重新整理,忽略空泡输运方程中的弛豫源项.整理后如下：

其中,U=(αρ,αρu,αρE)T,F=(αρ,αρu2+αP,αρE+αP)T,H=(0,pI,pIuI).

矿床内发育了两组断裂，一组为NNW向张扭性断裂，该组断裂以扭性为主，因扭力是顺时针方向的，故多为右斜列式出现，同时还伴生有小角度的扭裂隙和张裂隙。在这些构造裂隙带中充填有闪锌矿、辉锑矿、脆硫锑铅矿和石英、含锰菱铁矿、含铁菱锰矿、铁白云石等金属矿物和脉石。该组断裂多成群平行排列出现，地表出露30余条。它是本矿床的主要控矿构造，矿脉的产状及形态，严格受其控制。第二组为NNE向扭性断裂，分布在矿床西侧，多充填有乳白色石英脉，局部有辉锑矿化，有切穿NNW向断裂的现象，为第二期形成的断裂(图2)。在矿床主要地段未发现，对矿脉无破坏性。矿区内主要控矿构造及性质见表1。

求解格式如下：

式中：Δ=∂α1/∂x,其离散格式为

式中：S=max{|u i±c i|,|u i+1-c i+1|,|u i-1-c i-1|}.

因为空泡输运方程(16)为非守恒型方程,不能采用上述守恒形式格式求解.因此,对其采用如下求解格式：

1.2.2 边界条件

如图1所示,边界条件的构造采用特征线法(图中横坐标为网格长度上的网格节点,纵坐标为时间层).

图1 特征值及特征线在空间坐标的分布趋势d x/d t=u±c

相容性方程可以表示为：

式中：φR代表当特征波速为负值时,向上游传递的信息节点；φL代表当特征波速为正值时,向下游传递的信息节点；φA,φB,φC是不同网格节点上的值.

2 数值计算结果与试验结果对比

管线长度分为5个部分,管道长度在0～0.65 m之间,管道标高为0.47 m；管道长度在0.65～24.83 m 之间,标高为0.26 m；管道长度在24.83～49.17 m之间,标高为0.13 m；管道长度在在73.52～98.56 m之间,标高为0 m.计算时间步长为5×10-5s,管道波速为1 240～1 291.8 m/s.在管线沿线有4处压力检测点,分别位于管线1/4、2/4、3/4 和4/4 处.文献[12]给出了相应的流体热力学参数.使得该模型较好地保证了计算的稳定性和精度.

试验系统原理图如图2所示.

图2 试验系统示意图

表1 试验具体参数

数值模拟结果与试验结果对比如图3～8所示.

图3 管线P3处流速为2.74 m/s的水头随时间的变化特征

图4 管线P1处流速为2.74 m/s的水头随时间的变化特征

图5 管线P3处流速为1.43 m/s的水头随时间的变化特征

图6 管线P1处流速为1.43 m/s的水头随时间的变化特征

图7 管线P3处流速为0.47 m/s的水头随时间的变化特征

无论是管道末端还是管道中点的模拟结果都表明：采用非恒定摩阻的新模型,不仅可以捕捉到液柱分离过程中水击的最大升压和降压值,而且还可以捕捉最大升压之后的压力衰减变化.

3 内部流动特性分析

3.1 水锤二次升压分析

从图8可以看出,管道末端的水锤压力波的二次压力升高值明显高于首级波峰.这种现象是由液柱分离再弥合产生的压力波反射叠加所导致.

图8 管线P1处流速为0.47m/s的水头随时间的变化特征

为探究上述二次升压现象的流动机理,以“阀门-管道-稳压水箱”系统为研究对象,忽略摩阻影响.选取两种瞬变流动现象,第1种是初始流速为0.47 m/s,初始水头为70 m 无二次升压的水锤现象,其压力分布如图9(a)所示；第2种是初始流速为0.47 m/s,初始水头为50 m 有二次升压的水锤现象,其压力分布如图9(b)所示.

图9 阀门处两种瞬态流动压力分布

图10和图11为压力波在多相瞬变流中的传播过程,其中实线为升压波,虚线为降压波.压力波的传递过程遵循两个规律：(1)压力波与壁面、阀门等边界条件相接触时,其反射波与原波性质相同；(2)压力波与稳压边界相作用时,其反射波与原波性质相异.

图10 无二次升压的水锤现象

图11 伴有二次升压的水锤现象

如图10所示,压力波1、4、5为升压波,2、3、6为降压波.在无液柱分离情况下,压力变化按照理论计算公式(Hmax/min=H0±,H0为边界水头)下降与上升,压力下降与上升的变化范围如图12所示(横坐标0代表上游,1代表下游),该条件下并未出现二次升压现象.

图12 管线压力波分布

如图11所示,当末尾阀门快速关闭,阀门处产生升压波1向上游传播,升压波在上游与稳压罐作用形成降压波2反射至阀门处,传播时间为2L/c.阀门处的低压波2和3反射导致压力降低至饱和蒸汽压,阀门处发生空化.降压波3从上游稳压罐处变化为升压波4反射,升压波4在管道末端与空泡相遇,此时波速产生突变,如图13所示.

图13 波速随管线变化

升压波4在两相混合区域产生折射与反射两种压力波：第1种压力波与空泡作用反射形成降压波5向上游传递；第2种压力波在两相区域以升压波4继续向阀门边界折射传递.当压力波4传递至阀门处,反射为升压波6向上游传递.在降压波5和升压波6的传递期间,相同管线位置处降压波5 早于升压波6,其变化如图14(a)所示.

升压波6传递至上游处变化为降压波8向下游传递,降压波5传递至上游变化为升压波7向下游传递.在相同管线位置升压波7早于降压波8,其变化如图14(b)所示.升压波7在上游处与升压波6相遇,两种相同性质的正压波在N 点叠加传递至末端阀门处,最终使得末尾阀门处压力峰值高于理论计算值.

图14 管线压力波分布

3.2 水锤二次升压随系统参数变化

为了判定水锤二次升压值随系统参数的变化,本文选取了边界水头、边界流速、管道直径、管道长度、和水锤波速5个参数来进行内部流动特性分析.

3.2.1 边界水头对二次升压的影响

图15统计了二次升压值(二次峰值)在特定初始流速下随边界水头的变化,图中横坐标为边界水头(H0)与Joukowsky升压(ΔH=(cv0)/g)的比值,纵坐标表示二次升压值(H2max)与Joukowsky 升压值(ΔH)的比值.通过分析发现,在特定的初始流速下,二次升压值随着初始水头的增大呈现出先增大后减小的趋势,并且在H0≈0.8ΔH处会取得极值.此结论应用在水锤泵的输水过程中,可以使其获得更高的扬程.

图15 二次升压值随边界水头的变化

3.2.2 边界流速对二次升压的影响

表2统计了二次升压值在特定初始水头下随初始流速的变化.可以看出,在初始水头为24.8 m 条件下,初始流速在0.30～0.82 m/s之间,二次升压值随着初始流速的增大而升高,当初始流速大于1.02 m/s,二次升压现象便不再发生.在初始水头为41.8 m的条件下,初始流速为0.30 m/s,并没有发生二次升压现象,因为在此条件下的初始边界压力大于Joukowsky升压值,不能发生液柱分离,从而使得一次峰值与二次峰值相等.而初始流速在0.47～1.14 m/s之间,二次升压值随着初始流速的增大而增大,当初始流速大于1.35 m/s,二次升压现象便不再发生.在初始水头为71.8 m 条件下,初始流速在0.64～1.35 m/s之间,二次升压值随着初始流速的增大而增大,当初始流速大于1.70 m/s,二次升压现象便不再发生.以上说明,在液柱分离过程中,二次升压值随着初始流速的增大而增大,当初始流速大于临界流速,二次升压便不会发生.因此,在工程应用中,如果要避免因液柱分离而产生的二次升压,可以通过数值模拟来选择最优工况.

表2 二次升压值随初始流速变化(单位：Pa)

3.2.3 管道直径、管道长度和水锤波速对二次升压的影响

表3统计了在初始边界H0=0.8ΔH条件下,二次升压相对值随管道直径、管道长度和管道波速的变化.在理论情况下,压力波在水中的传播速度不超过1 450 m/s,实际波速值要小于理论值,并且为使得该条件下发生液柱分离,选取波速最小值为1 100 m/s.

表3 二次升压相对值H2max/ΔH 的变化情况

可以看出,随着管道直径的增加,二次升压相对值保持不变；随着管道长度的增加,二次升压值相应增大；随着水锤波速的增大,二次升压值随之增大.

3.2.4 二次升压试验验证

为了验证二次升压值随边界水头的变化规律,本节对图15的模拟结果进行了相应的实验验证.

表4 试验具体数据

首先进行了初始流速为0.348 3 m/s,边界水头为25.428 m 的水锤试验.在此基础上,保持初始流速基本不变,取边界水头为H0=0.6ΔH、H0=0.7ΔH、H0=0.8ΔH和H0=0.9ΔH,分别对相应试验2、3、4、5进行验证.试验结果如图16所示,可以看出,在H0=0.8ΔH处,二次升压值最大.

图16 二次升压值随边界水头变化

图17 管线中点处水头随时间的变化特征

上述试验结果对图15的模拟结果进行了验证.从水锤压力波随管线分布可以看出,由于摩阻的存在,升压波在阀门处反射向上游传递的过程中持续衰减,并且在管线中点处二次升压现象不再发生.

4 结论

本文通过对瞬态流动中的二次升压现象进行分析,得到了以下结论：

1)在伴有液柱分离的水锤中,空泡的生成与溃灭会使得管道内部发生压力波的叠加.管道中两种升压波叠加沿管线传递至阀门处,使得阀门处流体进一步被压缩,从而产生二次升压现象,这种二次升压值会高于理论计算值.

2)二次升压值在H0≈0.8ΔH处取得极值,随初始水头增大呈现先增大后减小的趋势.二次升压值随着初始流速的增大而增大,当初始流速大于临界流速,二次升压便不会再发生.二次升压值随水锤波速和管道长度的增大而增大,管道直径对二次升压值无影响.

3)通过5 组试验结果验证了二次升压值在H0≈0.8ΔH处取得极值.在升压波从阀门处向上游传递的过程中,管线阻力对其衰减影响显著.