在电厂、化工厂、化石燃料等各类输水管道系统中，由于水中气体释放、系统进气等原因，常会出现均匀分布的较小气泡(含气量小于1%)，从而形成水气两相均质流[1-4]．当阀门启闭或泵启停时，将会发生复杂的水气两相均质瞬变流现象，可能引发危险压力波动[5-6]，因此能精确模拟管道中均质流的瞬变现象具有重要意义[1]．文献[2]对水库-管道-水库系统中的两相均质流进行了试验研究，指出管道中的两相气液流动可以看作是小孔隙组分和均相混合物的等效单相流体，这也是后续均质流研究的基础；文献[3]通过比较数值结果和试验结果，评价了各种计算格式在气泡流和段塞流模拟中的准确性．有限体积法Godunov格式因能提高瞬态模拟的精度和效率近几年被广泛应用，文献[7]提出了一阶和二阶Godunov格式(GTS)来模拟两相瞬态流动，并证明了GTS方法在捕捉冲击时具有较高的分辨率．文献[8]针对两相水锤流提出了一种高效的GTS方法．与特征线、有限差分方法相比，二阶有限体积法Godunov格式具有计算精确、稳定、高效的优点，尤其是对于水气两相均质瞬变流模拟求解中常要处理库朗数小于1的情况更能体现该种方法的优越性．本研究采用有限体积二阶Godunov格式对水气两相均质流瞬变现象进行建模和模拟研究．为提高模型计算精度，考虑动态摩阻因素，将改进的Brunone动态摩阻模型引入数学模型中．采用二阶Godunov格式对数学模型进行求解；提出双虚拟单元的边界处理方法，实现了计算区域的所有节点和边界统一计算．将本研究所建数学模型计算结果与已有模型计算结果、试验数据进行对比，并对含气率、库朗数Cr等参数的敏感性进行分析．1 水气两相均质流模型1.1　基本控制方程水气两相均质流可视为具有平均性质的单当量流体[2]，其运动和连续性方程的矩阵形式为∂U/∂t+∂F/∂x=S，(1)U=ΩQm，F=QmQm/Ω+A1p，S=0(S0-S1)ρ1gA1，式中：ρ1为流体密度；A1为导管截面面积；Ω=ρ1A1为管道单位长度的流体质量；Qm=Ωv为质量流量，v为水流的流速；p为作用在截面A1重心上的压力；g为重力加速度；S0为导管的斜率；S1为摩阻项．式(1)中有Ω，p，Qm三个未知量，难以求解，引入压力波速的定义，有ag=[d(A1p)/dΩ]1/2．(2)气液混合的波速可表示为am=a/[1+ψrefρfrefa2pref1/β/p1+β/β]1/2，(3)式中：a为纯水体时水锤波速；pref为参考压力；ρfref为参考密度；β为多变指数，等温状态取1.0，绝热状态取1.4；Ψref为参考压强下含气率．在大气压下，4 ℃温度下压力、密度为参考压力和参考密度，取pref=10.325 kPa，ρfref=1 000 kg/m3．将式(3)代入式(2)并积分，得到Ω=Ωref+A1{p-pref+(pref-1/β-p-1/β)βψrefρfrefa2prref-1/β}/a2, (4)式中Ωref=ρfrefA1，p可以由式(4)迭代求出．若是纯水锤，则a恒定，保持不变．1.2　动态摩阻模型水力瞬变问题模拟计算常采用恒定摩阻近似计算瞬变流动过程中的摩阻项，虽然预测结果的最大值与试验符合较好，但无法准确模拟压力衰减[9-11]．为提高模拟精确度，考虑到改进的Brunone动态摩阻模型具有简单、计算快捷等优点，本研究将其引入控制方程．摩阻项分为稳态摩阻Js和非定常摩阻Ju两部分，即J=Js+Ju．Js通常根据水头损失的达西-维斯巴赫公式表示，即Js=fvv/(2gD)，其中：f为水头损失系数；D为管道直径．Brunone等人改进的动态摩阻模型，将非定常摩阻部分与局部瞬时加速度∂v/∂t和瞬时对流加速度∂v/∂x联系起来，即Ju=k(∂v/∂t+asign(v)∂v/∂x)/g，(5)式中：若v≥0，则sign(v)=1；若v0，则sign(v)=-1．式(5)中动态摩阻系数k=(C*)1/2/2，其中(C*)1/2为Vardy的剪切衰减系数，取决于流型．雷诺数Re=vD/υ，其中υ为运动黏滞系数．当Re2 320时，水流状态为层流，C*=0.004 76．当Re2 320时，水流状态为湍流，有C*=7.41/Relg(14.3/Re0.05)．2 二阶Godunov格式求解2.1　均质流的网格划分及黎曼求解器有限体积法[12]将空间域x离散为长度为Δx，时间域t离散为持续时间间隔Δt．第i个单元以节点i为中心，从i-1/2延伸到i+1/2．计算模型局部离散网格示意图如图1所示．10.13245/j.hust.211010.F001图1离散网格将式(1)沿x方向从控制边界i-1/2到控制边界i+1/2进行积分，并带入Ui=1Δx∫i-1/2i+1/2udx得到Uin+1=Uin+ΔtΔxFi-1/2n+1/2-Fi+1/2n+1/2+ΔtΔx∫i-1/2i+1/2Sdx，(6)式中：n表示t时刻；n+1表示t+Δt时刻；S为源项．本研究采用MUSCL-Hancock方法[12-13]进行二阶线性重构，从而获得空间为二阶精度的Godunov格式，为避免线性插值重构时出现虚假振荡的现象，引入斜率限制器函数MINMOD(σjn,σj-1n)[12]．对于水锤，流速v比a的数量级小的多，控制方程中可忽略对流项；然而am可能变为相对较小的值，须考虑对流项．当求解黎曼解时，以中间*区域的变量值U*来近似代替精确解．通过求解线性化双曲方程组的黎曼解得出Ω*=(ΩL+ΩR)/2+1+(vL+vR)/2；Qm*=QmL+(v¯-a¯m)(Ω*-ΩL)，式中：Ω*和Qm*分别为*区域的流体质量和流量；ΩL，QmL和vL分别为*区域左侧的流体质量、流量和流速；ΩR和vR分别为*区域右侧的流体质量和流速；v¯和a¯m分别为平均值．2.2　双虚拟单元的边界条件处理方法式(6)用于所有计算域的计算，包括上下游边界F1/2和FN+1/2，其中N为划分总单元数．每个Fi+1/2在计算时都须通过MUSCL-Hancock方法进行线性重构，此时须用到界面左右两侧的数据，但是划分控制体时只包含单元控制体1~N，因此文献[12]提出的二阶Godunov格式求解水锤问题在边界处实质是一阶精度．文献[7]提出在边界侧各加一个虚拟边界如图2(a)所示，边界处仍须单独处理．10.13245/j.hust.211010.F002图2计算区域与虚拟网格本研究提出了在边界侧增加两个虚拟单元，如图2(b)所示，使边界1和N能同时达到二阶精度，尤其是边界处与内部单元实现统一求解，更简单方便．2.2.1　二阶边界条件上一时刻计算域0~N+1内质量、流量均已知，并对单元1~N进行了二阶线性重构，根据黎曼不变量理论[5]，求得边界处的相关变量：(am/Ω)dΩ+dv=0    (沿着dx/dt=v+am)；(7)(am/Ω)dΩ-dv=0    (沿着dx/dt=v-am)．(8)以上游边界为例(下游边界类似)，沿dx/dt=v-am与单元1连接，对式(8)进行梯形积分近似得到(am1,Ln+ambn+1/2)(Ωbn+1/2-Ω1,Ln)/2-(Ωbn+1/2+Ω1,Ln)(vbn+1/2-v1,Ln)/2=0, (9)式中：ambn+1/2，Ωbn+1/2和vbn+1/2分别为左边界在n+1/2时刻的气液混合的波速、单位长度的流体质量及流速．am1,Ln+1/2，Ω1,Ln+1/2和v1,Ln+1/2分别表示单元1的左侧在n+1/2时刻的气液混合的波速、单位长度的流体质量以及流速．若已知上游边界水头hb，则压力pbn+1/2已知，Ωb可由式(4)求出；流速vbn+1/2由式(9)求得；vbn+1/2确定后，Qmbn+1/2也可求出．若已知上游边界流量Qb，则vbn+1/2已知；由式(9)可得Ωbn+1/2=1+2vbn+1/2-v1,Lnam1,Ln+ambn+1/2-vbn+1/2+v1,LnΩ1,Ln．(10)式(10)含Ωbn+1/2和ambn+1/2两个未知数，须通过以下步骤迭代求得：a. 假设Ωbn+1/2=Ω1,Ln；b. 将Ωbn+1/2代入式(4)求出pbn+1/2；c. 将步骤b求得的pbn+1/2代入式(3)求出ambn+1/2；d. 将步骤c求得的ambn+1/2代入式(10)求出Ωbn+1/2，重复上述步骤，求得Ωbn+1/2；e. 根据Qmbn+1/2=Ωbn+1/2vbn+1/2求出Qmbn+1/2．2.2.2　双虚拟单元的边界条件本研究在上、下游边界处外侧引入两个虚拟单元，其内参数与边界处参数相同，表示为Ω-1n=Ω0n=Ω1/2n=Ωbn；Qm-1n=Qm0n=Qm1/2n=Qmbn，式中：Ω-1n，Ω0n，Ω1/2n和Ωbn分别为单元编号-1，0，1/2和b处在第n时刻的流体质量；Qm-1n，Qm0n，Qm1/2n和Qmbn表示单元编号-1，0，1/2和b处在第n时刻的流量．2.3　动态摩阻项的时间积分与稳定性约束通量计算完成后，对式(6)进行积分，得到Uin+1=Uin+Δt(Fi-1/2-Fi+1/2)/Δx．采用二阶龙格-库塔进行时间离散，实现动态摩阻源项S的求解，得到如下显式过程：U¯in+1=Uin-Δt(Fi+1/2-Fi-1/2)/Δx；U¯'in+1=U¯in+1+S(U¯in+1)Δt/2；Uin+1=U¯in+1+ΔtS(U¯'in+1)．稳定性约束不仅包括对流部分的Courant-Friedrichs-Lewy(CFL库朗数约束)准则，而且包括源项的约束．由CFL得到ΔtCFL=min Δxi/(|vin|+|amin|)    (i=1,2,…,Nx)．二阶龙格-库塔离散的约束条件为ΔtS=min(-4Uin+1/S(Uin+1),-2Uin+1/S(U¯in+1))．综上所述，时间约束条件为Δtmax=min(ΔtCFL,ΔtS)．3 计算分析3.1　数学模型验证为验证所建数学模型，以文献[2]试验为依据，将本研究考虑动态摩阻的二阶Godunov格式水气两相均质流模型(2nd Godunov-UFM)模拟结果与相应恒定摩阻模型(2nd Godunov-SFM)的计算结果、文献[2]试验结果进行对比分析．案例1为文献[2]试验．文献[2]对管道内水气两相均质流动进行了细致的试验研究；试验装置如图3所示，管长为30.60 m，内径为0.026 m，初始上游水头Hr为21.7 m，流速为2.94 m/s，波速为715 m/s，水头损失系数为0.019 5，含气率为0.005 3，坡度为0．控制出口阀门和进口注入空气的压力，在管道中建立稳定的气-水混合物流动．当管道下游阀门快速关闭时，产生瞬态流动．管道中布有3个压力传感器，依次在距上游x1=8.0，21.1，30.6 m处．10.13245/j.hust.211010.F003图3试验装置图(m)试验过程中，保持上游水位恒定，下游边界处(x1=30.6 m)设置了压力传感器，将测得压力值作为下游边界值，如图4所示，图中：H为水头；t为时间．图5为压力传感器1和2处的计算和测量值，β=1.0，Cr=0.95，N=100．如图5所示，在上述试验条件下，本文数学模型能准确地模拟试验结果，包括压力峰值、周期和衰减．与仅加入恒定摩阻结果相比，加入动态摩阻计算结果更接近试验结果．10.13245/j.hust.211010.F004图4下游末端试验压力10.13245/j.hust.211010.F005图5x1=8.0，21.1 m处计算和试验结果图6为在等温(β=1.0)、中间过程(β=1.05)、绝热(β=1.4)条件下，本文的2nd Godunov-UFM模型计算值与试验值对比，Cr=0.95，N=100．与试验结果相比，绝热条件压力计算峰值偏大，这是因为试验中管道为不锈钢管，导电性高，并不满足绝热条件，瞬变过程中存在热交换；同时，等温条件压力计算峰值略小于试验值，这是由于瞬变时间较短(仅2.5 s)，很难进行充分的热交换而达到等温条件．所以，在文献[2]的试验中，绝对的等温和绝热条件都不存在．10.13245/j.hust.211010.F006图6x1=8.0，21.1 m处计算和试验结果实际上，文献[2]也指出：其试验热力学条件既不是等温的，也不是绝热的，是在等温和绝热条件1.0~1.4之间．经计算和试验对比发现：当β=1.05时，计算值与试验值整体符合最好，如图6所示．采用了相对误差平均值和均方根误差来定量评估模型精确性，如表1所示(案例1)．理论上，随着相对误差平均值和均方根误差值的减小，所提出的模型变得更加精确．表1的结果更加明确，可见：本文模型比只考虑恒定摩阻的模型更精确；当β=1.05时，计算值与试验数据整体符合最好．10.13245/j.hust.211010.T001表1误差分析表βx1/m模型相对误差平均值/%均方差1.008.0SFM5.771.90UFM5.481.8121.1SFM4.841.68UFM4.601.581.058.0SFM5.761.85UFM5.501.7821.1SFM5.311.62UFM3.901.551.408.0SFM9.092.49UFM11.272.9621.1SFM5.552.48UFM11.302.993.2　参数敏感性分析案例2为水库-管道-阀门系统：与文献[8]的案例相同，管长为10 km，内径为1 m，初始上游水头为200 m，流量为2 m3/s，参考压力为101.3 kPa，含气率为0.2%；初始波速为1 000 m/s．下游阀门瞬间关闭，为避免虚假数值震荡，库朗数取为0.95，N=200，β=1.0．文献[8]假设管道光滑无摩阻，本研究在此基础上考虑管道存在摩阻，并加入动态摩阻．图7给出当t=140 s，网格数为200时，管道沿线压力计算结果．图7中x2为纵向距离．图8给出了管线中点处瞬变压力计算结果．图7和8均显示：与本文模型计算结果相比，无摩阻模型所计算的压力值及发生时间均存在较大误差；仅考虑恒定摩阻模型所计算的压力值略高，但时间相位差较大．10.13245/j.hust.211010.F007图7案例2纵向压力变化图10.13245/j.hust.211010.F008图8案例2管线中点处瞬变压力a.　含气率的影响图9显示了案例2中初始含气率分别为0.01，0.005 3，0.002，0.001和0时管线中点处瞬态压力．结果表明：随着含气率的增加，波速减小，压力峰值减小，振荡周期增大；同时，当含气率为0时，是纯水锤问题，用文献[12]发表的水锤的二阶Godunov求解格式加入Brunone动态摩阻(边界处理用本文的双虚拟单元的边界处理法)同本文含气率为0的水气两相均质流模型进行对比，发现数据重合度高，再次证明了模型的正确性．10.13245/j.hust.211010.F009图9含气率对管道中间截面水头的影响b.　Cr的影响采用案例2研究不同Cr的影响，图10分别给出了管线中点处、无摩阻和考虑动态摩阻的压力计算结果．图10(a)的结果表明：当Cr=1.0时，在最低压力附近，随着压力的急剧上升，出现了轻微的虚假数值振荡，这是由于当Cr接近或等于1.0时，在瞬变过程中压力突然增大，波速急剧上升，导致Cr可能大于1．因此，当Cr接近或等于1.0时，在陡坡附近会产生初始数值不稳定．当Cr=0.95，0.50，0.10时，可以消除激波引起的数值不稳定问题．图10(b)结果显示：本文数学模型在Cr接近或等于1.0处仍能保证计算结果的稳定性和准确性，可有效避免虚假的数值震荡．10.13245/j.hust.211010.F010图10Cr对瞬变压力影响总体上，较低的Cr会导致一定的数值耗散；同时，随着Cr的降低，计算效率迅速降低，因为较小Cr导致较小的时间步长．4 结论a. 与现有的均质流模型相比，本研究建立的考虑动态摩阻二阶Godunov格式水气两相均质流模型能更准确地模拟试验压力波动曲线．b. 水气两相均质流快速瞬变过程中，微量气体压缩膨胀过程并非绝热过程而接近等温过程，建议多变指数取值1.05．当多变指数取值1.4时(绝热)，瞬变压力计算结果与实测值存在较大差别；当多变指数取值1时(等温)，瞬变压力计算结果与实测值基本符合；当多变指数取值1.05时，瞬变压力计算结果与实测值更为接近．c. 含气率对瞬变压力具有重要影响．随着初始气体含气率增加，水锤波速降低，瞬变压力峰值减小，压力周期延长，因此输水管道中存在一定量自由气体可缓解水锤压力．d. 动态摩阻的引进提高了计算结果稳定性，保证了计算精确度．当库朗数接近1时，不考虑动态摩阻的数学模型可能会出现数值振荡的计算不稳定现象；考虑动态摩阻的数学模型能有效避免上述虚假数值振荡，提高计算结果稳定性和精确度．