在长距离调水工程、水电站、泵站、火电核电冷源系统等有压输水系统中，系统启动、停运、调控等因流量变化而发生水力瞬变(或水锤)，若水力元件操作或水锤防护措施不合理，则会引起异常的压力波动，导致管线结构破坏甚至爆管 [1-4]．因此，准确模拟压力管道中的快速瞬变，对输水系统安全稳定运行和实现智慧精准化调度运行极为重要．目前，在有压输水系统水力过渡过程计算中，通常考虑稳态摩阻或准稳态摩阻的影响，能快速地模拟水力瞬变过程，预测出水锤发生过程中的极值．但实际水锤发生过程中能量耗散受到各种复杂因素影响，其压力计算结果与实测数据存在一定偏差，尤其是低估了瞬变压力的衰减[5]．实际输水管线长、特性变化多，水锤波相互传递叠加的过程错综复杂，准确模拟水锤波的衰减过程对于后续水力安全评估、实时监测分析非常重要．为了准确模拟瞬变压力衰减效应，须要将动态摩阻引入水锤控制方程中，主要有瞬时加速度动态摩阻模型(即Brunone模型[6])和加权类动态摩阻模型．Brunone类模型构成简单，其易于各类求解格式融合，但属于经验性模型．加权类模型被公认为计算相关摩阻损失的精确模型．Zielke[7]基于层流推导水锤发生过程中的动态摩阻模型，其能准确模拟管道层流和低雷诺数紊流的瞬变过程．文献[8-10]推导出适用于紊流的加权类动态摩阻模型．为了提高加权类动态摩阻模型计算效率，大量研究者对模型修改进行了研究，提出了高效简化卷积模型，例如Trikha-Vardy-Brown (TVB)[11]模型．虽然动态摩阻模型能提高水锤模拟精度，且计算效率得到了改进和提升，但在实际输水系统中依旧很少应用；缺乏动态摩阻对水锤压力衰减效应的实验验证数据，动态摩阻衰减效应及运行参数影响规律亟待深入研究，以便准确地模拟瞬变过程．本研究设计搭建了水锤实验系统，验证了高效简化卷积动态摩阻模型TVB的准确性；对比分析稳态、动态摩阻模型模拟水锤的波峰差异大小规律，提出了描述动态摩阻相对重要性的无量纲参数K，分析了管道长度直径比、马赫数、管道摩擦系数、初始雷诺数等敏感性参数对水锤压力衰减的影响．1 含动态摩阻的水锤模型及其求解1.1　水锤控制方程对于管道内的有压流，其控制方程可由连续性方程与动量方程得到，可写成如下的经典水锤方程∂H∂t+a2g∂V∂x=0；(1)∂H∂x+1g∂V∂t+J=0，(2)式中：H为测压管水头；V为流速；g为重力加速度；a为波速；J为管道流体中由于摩擦引起的单位长度的水头损失；x为沿管轴线距离；t为时间．1.2　原始Zielke模型Zielke认为管道瞬变流中瞬时管壁切应力由恒定项和附加项组成，恒定项为恒定状态切应力值，产生的压力衰减与稳态摩阻相同．附加项利用加权函数考虑了历时速度和加速度对当前流态的影响，故摩阻项由恒定和非恒定摩阻项组成，有J=JS+JU，(3)式中JS和JU分别为管道流中由恒定摩阻和非恒定摩阻引起的水头损失．恒定摩阻项JS=fVV/(2gD);(4)非恒定摩阻项JU=16νgD2∫0t∂V∂t(u)W(t-u)du,(5)式中：ν为运动黏滞系数；u为在0～t之间积分的时间变量；f为达西威斯巴赫管道摩擦系数；D为管道直径；W为关于无量纲时间τ的加权函数，τ=4νt/D2．W具体计算如下．当τ≥0.02时，W(τ)=∑i=15exp(-miτ)；(6)当τ0.02时，W(τ)=∑i=16niτ(i-2)/2，(7)式中：mi，ni为系数，mi={26.374 4，70.849 3，135.019 8，218.921 6，322.554 4}，ni ={0.282 095，-1.250 000，1.057 855，0.937 500，0.396 696，-0.351 563}．1.3　Vardy-Brown模型Zielke模型能够准确模拟管道层流和低雷诺数紊流的瞬变过程，对于高雷诺数的紊流流动，Vardy和Brown提出了与雷诺数和时间相关的紊流加权函数，满足2 000Re1×108时的紊流模拟．改进的加权函数为W(t-t')=αexp[-β(t-t')]/π(t-t'),(8)式中：α=D/(4ν)；β=(0.54νkR)∕D2，kR=Relg(14.3∕Re0.05)，Re为雷诺数．1.4　高效简化卷积模型原始Zielke模型被公认为是准确的，但计算所需储存空间较大，执行起来很耗时．为了提高计算效率，很多研究者采用新的函数对加权函数进行近似代替，对Zielke模型进行了简化．Trikha首次对Zielke模型加权函数进行简化，提出一种近似加权函数，能够准确地模拟管道层流水力瞬变过程中的水锤压力、压力波变形及压力衰减，并且较大程度地缩短了计算时长，提高了模型的计算效率．Vardy和Brown拓展了Trikha提出的近似加权函数，在紊流流态下也能准确高效模拟管道中的瞬变过程．拓展的近似加权函数形式为Wapp(τ)=∑i=117mie-niτ，(9)式中mi，ni为系数，mi={9.06，-4.05，12，8.05，22.7，35.2，65.9，115，206，365，651，1 150，2 060，3 630，6 640，10 700，26 200}，ni={10，1×101.5，1×102，1×102.5，1×103，1×103.5，1×104，1×104.5，1×105，1×105.5，1×106，1×106.5，1×107，1×107.5，1×108，1×108.5，1×109}．本研究采用TVB高效简化卷积模型，使用矩形特征线网格代替Zielke的交错网格，通过卷积积分的一阶近似来计算动态摩阻，在保证水锤计算准确性的基础上，同时满足计算效率的要求．有    JU=16νgD2∑j=1k-1Vi,j+1-Vi,jWk-jΔt-Δt2=16νgD2∑j=1k-1Vi,k-j+1-Vi,k-jWjΔt-Δt2.(10)2 实验装置及工况为了验证本文模型在不同规模管道系统中模拟水力瞬变的正确性，在两个实验系统进行实验验证(见图1和图2)．两个实验系统均由上游的压力罐、潜水泵、电磁流量计、手动阀门、压力传感器、中间段以环绕形式固定的供水管道及下游水箱组成．实验系统一的管道是一根长241.52 m的钢管，管道内径为50 mm．实验系统二的管道是一根管长10.13245/j.hust.240285.F001图1水锤实验系统一10.13245/j.hust.240285.F002图2水锤实验系统二1 149.3 m的铜管，管道内径为20 mm，管道中间设置一个控制阀门．在实验中通过控制阀门的快速关闭产生水锤现象，阀门处设置的压力传感器记录水力瞬变时的压力变化．实验采用的高精度压力传感器BYP300的主要技术参数为：输出为4～20 mA，精度为0.2%，测量范围为0.0～1.0 MPa．电磁流量计的测量范围为0.707～10.000 m3/h．同时使用采样频率为2 kHz的LabVIEW软件进行数据采集．实验涉及到的尺度范围较广，其中管道长度直径比为4.83×103～5.74×104，雷诺数为2.16×103～1.77×104，包含层流和紊流两种流态．实验工况基本参数见表1，阀门处压力传感器测得的压力变化结果与数值模拟结果见图3．其中实验管道均为光滑管道，水体运动黏度取1×10-6 m2/s．10.13245/j.hust.240285.T001表1实验工况基本参数工况L/mD/mfa/(m∙s-2)V/(m∙s-1)1241.520.050.0141 3280.2802241.520.050.0141 3280.3543574.640.020.0301 2800.3014574.640.020.0301 2800.10851 149.300.020.0301 2800.30561 149.300.020.0301 2800.11610.13245/j.hust.240285.F003图3关阀水锤过程中稳态、动态摩阻模型压力计算结果与实验结果对比3 模型验证如图3所示，针对水锤实验系统一(241.52 m的钢管)的工况1和工况2，稳态摩阻模型能准确计算第一峰值水锤压力，但是从第二个波动周期开始，虽然压力波动周期与实测值基本一致，但其模拟出的压力峰值远高于实测值，偏差随时间的增加而不断增大．针对水锤实验系统二(574.64，1 149.30 m的铜管)的工况3～6，稳态摩阻模型计算出的第一峰值水锤压力略低于实验值，其原因是稳态摩阻模型没有考虑动态剪切应力在管壁水体流向反转时对水锤压力的正向作用[12]．从第二个波动周期开始，虽然压力波动周期与实测值基本一致，但其模拟出的压力与实测值存在较大差异，偏差随时间的增加而不断增大．在不同管道规模的水力瞬变实验中，针对6种水锤实验工况，本研究采用的TVB动态摩阻模型均能准确模拟出各个波动周期的压力峰值、波动周期、压力衰减量，与实验结果保持一致，且在整个计算时间域内均能较好反映水锤波的瞬变过程．此外，对比6个实验工况中稳态、动态摩阻模型模拟水锤波的差异大小，即动态摩阻的影响效应，与管道规模、初始流速等参数密切相关．表2给出了6种实验工况下稳态、动态摩阻模型两种模型压力计算误差与计算耗时对比，表中Pi为第i个水锤周期压力峰值的计算误差，Ti为第i个压力峰值出现时间的计算误差．结果表明：TVB模型能准确模拟不同管道规模的瞬态压力，前三个压力峰值最大误差不超过3%，验证了本研究建立的TVB模型计算准确性．此外，尽管TVB模型提高了加权类动态摩阻模型的计算效率，但与稳态模型相比，计算效率依旧低下(6种实验工况两者计算时间最小相差11倍)．对于压力峰值出现时间的计算误差，两种模型之间的差异并不是很明显(误差均小于2%)．10.13245/j.hust.240285.T002表2不同实验工况下稳态、动态摩阻模型压力计算误差与计算耗时参数摩阻模型实验工况123456Tcpu/s稳态0.3750.2810.2660.2660.2190.188TVB4.2664.3594.1414.0633.0312.813P1/%稳态0.0000.4681.8491.4992.9691.676TVB0.0000.4680.6870.0000.7940.386P2/%稳态2.2531.9320.2480.7241.8181.089TVB1.2480.0000.0030.0000.9090.508P3/%稳态4.1573.9802.6662.7782.5422.590TVB0.2310.4970.2660.0002.5421.004T2/%稳态0.1290.0001.2970.4811.3050.791TVB0.0000.0000.1180.0000.8460.000T3/%稳态0.1830.6230.4670.8650.2901.549TVB0.1570.0000.1480.0671.9500.3324 动态摩阻压力衰减理论分析4.1　无量纲分析非恒定项是加权类动态摩阻在水锤波衰减方面与稳态摩阻产生差异的根本原因，为了研究系统初始参数(管径和管长等)对动态摩阻非恒定项的影响，在水锤基本控制方程中推导初始参数对瞬态压力衰减的影响是有必要的．对系统水力参数进行无量纲化如下    x*=x/L;V*=V/V0;t*=t/(L/a);H*=H/H0;α*=α/D2/ν;β*=βD2/ν;α'=α*aD2/(νL);β'=β*[νL/(aD2)];W(t*)=α'e-β't*/πt*, (11)式中：L为管道长度；V0为管道初始流速；H0为上游水库水位．将无量纲参数代入到水锤方程(1)和(2)进行无量纲化，可得：∂H*∂t*+aAV0gH0∂V*∂x*=0；(12)    ∂V*∂t*+gH0aAV0∂H*∂x*+fLV0aDV*+4πvLa∫0tW(t'*)∂V*∂t*dt'*=0.(13)利用傅里叶变换，方程(12)和(13)可变为：iw*H*+aAV0gH0∂V*∂x*=0；(14)    iw*V*+gH0aAV0∂H*∂x*+fLV0aDV*+4πvLaiw*α'β'+iw*V*=0,(15)式中V*(x*,w*)和H*(x*,w*)为频域中管道流速和水头的振幅大小．联立式(14)和(15)，消除V*可得：∂2H*/∂(x*)2+CHH*=0；(16)    CH=(ω*)2[1-i[fLV0/(aDω*)]+[4πνkL/(aA)](α'/β'+iω*)].(17)上述方程的解形式为H*(x*,ω*)=H0*exp(-Kx*ω*)，(18)式中：H0*为频域中水头的振幅；K=Kr+iKi，Kr和Ki分别为水锤波衰减和相位变化的系数．4.2　描述动态摩擦相对重要性的无量纲参数Ghidaoui等[13]提出无量纲参数P，其物理意义是水锤波的径向扩散时间与轴向时间尺度的比值，P=2D/(Vf)/(L/a)．(19)参数P能评估水锤水力瞬变准稳态假设的准确性，也可以描述水锤波的衰减程度，但没有从理论上进行解析推导．文献[14]对方程(17)进一步简化后推导了水锤波衰减系数的表达式，其中Krs和Kru分别为恒定摩阻项和非恒定摩阻项的衰减大小，有    Kr=Krs+Kru=(fMa/2)(L/D)+2(L/D)(Ma/Re0)    (Ma(L/D)≪1);(fMa/2)(L/D)+[Re0/(0.54kR)1.5]∙[D/(LMa)]    (其他),(20)式中Ma为马赫数，Ma=V/a．Duan用无量纲参数I来描述动态摩擦的相对重要性，Wahba[15]利用类似形式的参数λ(λ=0.5I)研究水锤瞬态衰减程度，这些研究指出参数I可以描述动态摩擦的相对重要性，I=fMa(L/D)．(21)Meniconi等[16]分析证实，初始雷诺数Re0对水锤波的衰减同样有显著影响．管道系统的初始雷诺数Re0越大，稳态摩阻模拟的压力曲线与实验结果越符合，侧面反映了动态摩阻的相对重要性随初始雷诺数的增加而降低，说明Re0同样是描述动态摩擦相对重要性的指标．在有压输水系统中，为了更好地描述管道系统水力瞬变过程中动态摩阻的相对重要性，本研究提出新的无量纲参数K来描述动态摩擦相对重要性，K=k(Re0)/I.(22)k是只关于初始雷诺数Re0的函数，在后文研究初始雷诺数Re0对动态摩阻相对重要性的影响中，对结果进行拟合，得到k=Re0-0.45．(23)无量纲参数K可表示为K=Re0-0.45D/(fMaL)．(24)4.3　评价标准为了直观地衡量水力瞬变过程中动态摩阻对水锤衰减的影响，本研究采用局部瞬态分析的方法，采用参数γ来衡量动态摩阻非恒定项对水锤峰值的影响，从而评价动态摩阻对水力瞬态过程的相对重要性，γ(z)=[ΔHN(z)-ΔHT(z)]/ΔHN(z)，(25)式中：γ为非恒定摩阻项在动态摩阻水头损失的占比，γ越小表示动态摩阻非恒定项影响越小，动态摩阻与稳态摩阻模拟结果差异越小；ΔH为水锤波峰压力相对于初始波峰的压力差(即初始峰值压力减去后续峰值压力)，下标N表示TVB模型模拟结果，下标T表示稳态摩阻模型模拟结果；z为压力波峰数．5 参数敏感性分析通过动态摩阻压力衰减理论分析可知，动态摩阻的相对重要性取决于管道系统的长度直径比L/D、马赫数Ma、管道摩擦系数f和初始雷诺数Re0．为了分析这些参数的改变对瞬变压力峰值的影响，对这些参数相对变幅一致的不同取值进行模拟．TVB模型的准确性通过实验得到了验证，在此以TVB模型结果为参照，通过对比稳态摩阻模型模拟结果，分析动态摩阻模型的敏感性参数，探究动态摩阻相对重要性随敏感性参数变化的规律．本研究设置的基本计算工况如图4所示，在水库-管道-阀门系统中，上游水头Hr保持100 m，管长L=5 000 m，管径D和管道摩擦系数f按实际输水管分别取为1 m和0.015，初始流速V0＝0.5 m/s，阀门瞬间关闭以形成瞬变流动．计算网格数量设置为150，满足库朗数条件并有足够的精度．整体以控制变量的形式设置不同工况，从峰值衰减情况对摩阻模型准确度进行分析，并进一步考察关键参数对两类摩阻模型的模拟影响程度．10.13245/j.hust.240285.F004图4水库-管道-阀门系统图5.1　长度直径比的影响在有压输水工程中，输水管道的长度直径比L/D一般在1×102～1×105的尺度区间，为了衡量不同规模的输水系统分别采用两种摩阻模型对水力瞬变过程的模拟效果，通过改变管道长度来控制L/D的取值，本研究分别选择L/D=5×102，1×103，5×103，1×104，5×104，1×105的工况进行模拟，并绘制前50个水锤波周期的γ变化如图5所示．10.13245/j.hust.240285.F005图5前50个周期γ- L/D的变化图从图5可知：随着L/D的增加，动态摩阻非恒定项造成的压力衰减越来越小，动态摩阻的相对重要性逐渐降低．利用前10个周期平均γ进行分析(不含首个周期)，当L/D=500时，前10个周期γavg=0.861，稳态摩阻模型产生的压力衰减仅为动态摩阻模型的13.9%，此时不可忽略动态摩阻对压力衰减的贡献．当L/D=5×104时，γavg=0.008，稳态摩阻模型产生的压力衰减为动态摩阻模型的99.2%，稳态与动态摩阻在压力衰减模拟精度上基本一致．同时，由图中γ的变化可知，随着水锤波在管道中不断传播，非恒定摩阻项造成的压力衰减程度是逐渐减小的．将不同L/D对应的γavg结果进行拟合，得到其拟合规律为γavg=4 908/(L/D+4 806).(26)由拟合公式可知，在其他敏感性参数确定的情况下，γavg与L/D存在某种关系，由拟合结果可认为动态摩阻的相对重要性与L/D大致呈反比例关系．5.2　马赫数的影响在输水管道系统中，马赫数的取值范围一般是1×10-4～1×10-3的尺度区间，为了研究马赫数对动态摩阻模型在水力瞬变中相对重要性的影响，设置了不同波速大小从而控制马赫数取值，分别选择Ma=1.25×10-4，2.50×10-4，5.00×10-4，1.0×10-3，2.5×10-3的工况进行模拟，并绘制前50个周期γ的结果如图6所示．10.13245/j.hust.240285.F006图6前50个周期γ-Ma的变化图从图6可知：随着Ma的增加，动态摩阻非恒定项造成的压力衰减越来越小，动态摩阻的相对重要性逐渐降低．当Ma=1.25×10-4时，γavg=0.836，稳态摩阻在前10个周期(不含首个)产生的压力衰减仅为动态摩阻的16.4%．当Ma=2.50×10-3时，γavg=0.071，压力衰减几乎都由恒定项产生，稳态摩阻模型产生的压力衰减占动态摩阻模型的92.9%．将不同M对应的γavg结果进行拟合，得到其拟合规律为γavg=4.812×10-4/(Ma+4.144×10-4)．(27)由拟合公式可知，在其他敏感性参数确定的情况下，γavg与Ma存在某种关系，由拟合结果可认为动态摩阻相对重要性与Ma大致呈反比例关系．5.3　管道摩擦系数的影响为了研究管道摩擦系数f对动态摩阻模型在水力瞬变中相对重要性的影响，分别选择管道摩擦系数f=0.005，0.010，0.015，0.030的工况进行模拟，并绘制前50个周期γ的结果如图7所示．10.13245/j.hust.240285.F007图7前50个周期γ-f变化图从图7中可以看出：随着管道摩擦系数f的增加，动态摩阻非恒定项造成的压力衰减越来越小．当f=0.005时，γavg=0.81，稳态摩阻在前10个周期(不含首个)产生的压力衰减仅为动态摩阻的19%；当f=0.030时，γavg=0.38，稳态摩阻产生的压力衰减为动态摩阻的62%，动态摩阻在水力瞬变中相对重要性随f的增大而减小．将不同f对应的γavg结果进行拟合，得到其拟合规律为γavg=0.022 49/(f+0.002 263)．(28)由拟合公式可知：在其他影响参数确定的情况下，γavg与管道摩擦系数f存在某种关系，由拟合结果可认为动态摩阻重要性与f大致呈反比例关系．5.4　初始雷诺数的影响为了研究初始雷诺数Re0对动态摩阻模型在水力瞬变中相对重要性的影响，通过改变D且保持L/D不变，从而控制初始雷诺数的取值，分别选择初始雷诺数Re0=2.50×105，5.00×105，1.00×106，1.67×106的工况进行模拟，并绘制水锤波前50个周期的γ结果如图8所示．10.13245/j.hust.240285.F008图8前50个周期γ-Re0变化图从图8中可以看出：在水锤波的前50个周期，随着初始雷诺数Re0的增大，动态摩阻非恒定项造成的压力衰减并不是越来越小，原因主要是由于动态摩阻模型采用的基于Re0条件的冻结紊流假设在瞬态早期阶段最有效[13]，随着时间的推移，动态摩阻模型模拟的水锤波开始失去准确性．对于阀门突然关闭的瞬态，管道流体的平均流速和紊流结构会随时间发生改变，使得紊流的黏度分布规律和剪切层厚度发生变化，流动将逐渐失去对Re0的依赖性，导致计算精度随着模拟时间的增加而降低．考虑水锤波的前10个周期，Re0造成水锤波衰减的规律相同，利用γavg来分析Re0对水锤压力衰减的影响(见图9)．当Re0=2.5×105时，γavg=0.614，稳态摩阻产生的压力衰减为动态摩阻模型的38.6%；当Re0=1.67×106时，γavg=0.426，稳态摩阻产生的压力衰减为动态摩阻模型的57.4%．随着初始雷诺数Re0的增加，动态摩阻非恒定项造成的压力衰减程度越来越小．将不同Re0下的γavg结果进行拟合，得到其拟合结果为γavg=235.9Re0-0.45．(29)10.13245/j.hust.240285.F009图9前10个周期γ-Re0变化图由拟合公式可知：在其他影响参数确定的情况下，γavg与Re0存在某种确定关系，由拟合结果可认为动态摩阻重要性与Re0大致呈幂函数关系．5.5　无量纲参数K的分析影响动态摩阻相对重要性的敏感性参数与压力衰减理论分析结论一致．计算表明以下无量纲参数K可用于描述有压水力瞬态中动态摩阻的相对重要性，K=Re0-0.45D/(fMaL)．(30)在实际有压输水工程中，无量纲参数K的取值范围一般不超过1，不同工况下的γavg随无量纲参数K变化的结果如图10所示，由图10可知无量纲参数K与γavg的大小保持高度一致性，γavg随无量纲参数K正向变化．无量纲参数K可用于描述有压水力瞬态中动态摩阻的相对重要性，K值越大，动态摩阻在有压水力瞬变压力衰减中越重要．10.13245/j.hust.240285.F010图10γavg随无量纲参数K变化图由图10可知：当0.0 K ≤ 0.1时，γavg对K值大小特别敏感，此时动态摩阻对压力衰减的贡献随K的改变而显著变化；当0.1 K ≤ 1.0时，γavg的大小保持一个较高水平，此时不可忽略动态摩阻在水力瞬变压力衰减中的贡献．利用曲线拟合，图10的结果可以近似为γavg=70K2+14.82K-0.092 1(0.0K≤0.1);-0.014 9K-1.166+0.893 1(0.1K≤1.0).6 结论a．本研究建立了考虑动态摩阻效应的TVB水锤模型，同时设计搭建了两个水锤实验系统，实验验证了TVB动态摩阻模型的准确性．在不同管道规模及广泛雷诺数的水力瞬变中，TVB模型均能准确地描述压力波的波形畸变和衰减过程．b．通过对动态摩阻压力衰减的理论分析，确定了动态摩阻相对重要性的敏感性参数，揭示了动态摩阻的影响效应：管道直径长度比、马赫数、管道摩擦系数越大，初始雷诺数越小，动态摩阻在水力瞬变过程中对水锤压力衰减的相对贡献越小．c．提出了能准确描述动态摩阻相对重要性的无量纲参数K，K值越大，动态摩阻在有压水力瞬变过程中对水锤压力衰减的相对贡献也越大．该无量纲参数可为有压输水系统的设计优化及安全运行提供支撑．