水电能源是可再生能源的最大来源[1]，中国是目前世界上最大的水电生产国[2]．在高水头多泥沙流域，当水头达到300~400 m以上时，水斗式水轮机在寿命周期成本上更加具有优势[3-4]．转轮与配水环管是水斗式水轮机的核心部件，葛新峰等[5-8]研究了泥沙磨损及磨蚀对转轮的影响．配水环管是水斗式水轮机流量分配及导流的重要部件，位于进水阀和喷嘴之间，可以有效引导水流使其呈轴对称进入各个岔管及其喷嘴[9]．针对国内某水电站在电网中的重要性，对各方案各种可能工况下的过渡过程进行数值模拟计算和分析，并根据计算结果和已有的设计理论、设计经验对输水系统进行优化，以便为输水系统布置、结构设计等工作提供可靠依据．配水环管的内部流动特性及压力变化直接影响到机组运行安全，特别是在甩负荷过程中其瞬态流动特性比较复杂．对于配水环管的流动问题，大多数都是依托模拟软件进行压力和流动研究．在流动方面，陈创新等[10]研究了不同开启状态下的配水环管中的流动规律及流动损失，Han等[11]研究了二相流中的配水环管及分岔管的流动现象；在压力方面，樊小发等[12]利用充水保压技术对配水环管的安装提供了一定的理论支撑．对于过渡过程问题，李金伟等[13]对混流式水轮机进行飞逸过渡过程的三维非定常湍流计算，刘梦瑶等[14]利用DEFINE_CG_MOTION宏进行简单的网格运动定义．本研究建立国内某冲击式电站的配水环管、喷针及折向器模型，采用水-气两相模型进行了非定常计算，利用动网格方法研究过渡过程中六喷嘴配水环管产生的流速及压力的动态变化．1 计算模型国内某电站额定水头671 m，引水线路长5.5 km，配水环管采用6喷嘴的月牙肋岔管设计，进水口直径为2 700 mm，喷嘴口直径为400 mm，配水环管出口直径为1 250 mm，喷针行程为370 mm，喷针关闭规律为直线关闭25 s，初始开度T1=0.75．配水环管计算域的三维几何模型及其截面、引水线路布置如图1所示．10.13245/j.hust.240653.F001图1基础模型及调速机构关闭规律喷针开度T1通常设定为剩余行程与总行程的比值，决定了喷嘴出口的流量．折向器偏转角度为30°，动作时间为3 s，关闭规律为直线关闭，不同折向器开度(T2)示意图如图2所示．10.13245/j.hust.240653.F002图2不同折向器开度示意图2 数学方程及计算方法2.1　数学模型2.1.1　湍流模型在模拟喷针关闭过程中，内部复杂流动时冲击式水轮机内部流动为弯管湍流，因此采用标准k-ε模型，其稳定性和计算精度均比较高，适用于模拟充分发展的湍流，能较好地预测配水环管内部复杂流动特征[15]，湍动能方程和湍动能耗散率方程如下．湍动能方程(k方程)∂(ρk)∂t+∂(ρkui)∂xi=∂∂xjμ+μtσk∂k∂xj+Gk+Gb-ρε-YM+SK，湍动能耗散率方程(ε方程)∂(ρε)∂t+∂(ρεui)∂xi=∂∂xjμ+μtσs∂ε∂xj+G1εε(Gk+G3εGb)/k-G2ερε2/k+Sε，式中：xi和xj为i和j点坐标位置；ui为i方向的速度；σk为k方程的湍流普朗克数；t为时间；ε为湍动能耗散率；σs为ε方程的湍流普朗克数；μt为湍流黏性系数，μt=ρCμk2/ε(ρ为密度；Cμ为模型常数，Cμ=0.09；k为湍动能)；Gb为由浮力引起的湍动能生成项；Gk为由平均速度梯度引起的湍动能生成项；YM为由可圧缩湍流中的脉动产生的湍动能生成项；C1ε，G2ε，G3ε和G4ε为经验常数，C1ε=1.44，G2ε=1.92，G3ε=0.99，G4ε=0.09；SK和Sε为用户定义的源项．2.1.2　多相流VOF模型采用的多相流模型为VOF(volume of fluid)气液两相流，其控制方程如下1ρq∂∂t(aqρq)+∇⋅(aqρqv⃗q)=0；∑q=1naq=1；aqn+1ρqn+1-aqnρqnΔtV+∑fρqn+1Ufn+1aq.fn+1=0；ρ=∑q=1naqρq；∂∂tρv⃗+∇⋅(ρv⃗v⃗)=-∇p+∇μ∇v⃗+∇v⃗r+ρg+F，式中：n为上一个时间步的索引；Δt为计算步长；Ufn+1为第n+1步时通过表面的体积通量；aq.fn+1为第n+1步时体积分数的面值；F为表面张力的等价体积形式；∇为哈密顿算子；g⃗为重力加速度；ρq为第q相的物理密度；∇v⃗r为平均速度矢量；v⃗q为第q相的速度；aq为第q相的体积分数；ρ为密度．2.2　计算域和网格划分及有效性验证喷嘴末端延伸出一段区域作为射流的空气域，针对瞬态模拟，本次计算采用压力进出口条件．按照电站水流在配水环管进口处的总压力水头为580 m，换算成进口压力为5.69 MPa，湍流参数I=2%，水力直径为1.721 m．配水环管作为复杂几何外形的管道设备，采用四面体单元网格进行非结构网格的划分方法[10]．由于喷针及折向器在计算过程中处于运动状态，会极大地对网格进行拉扯变形，甚至在阀门即将关闭时，网格尺寸会无限趋于0，因此必须对喷嘴喉部及周边的计算区域进行网格加密处理，防止因为网格变形幅度过大导致计算失败．为了检验网格的独立性，采用了5种不同的网格数(n=1.537 755×106，2.733 022×106，3.467 954×106，3.869 673×106，4.434 595×106)，网格划分如图3所示．10.13245/j.hust.240653.F003图3配水环管网格划分对划分的5种不同方案，在数值模拟计算得到射流中心的速度值之后，与计算理论值进行比较．当n=3 869 673时已经能满足需求，增加网格数会导致计算时间的大幅度上升，因此选用总网格数为3 869 673．2.3　动网格及有效性验证2.3.1　动网格算法与UDF代码编制喷针的关闭过程属于流体工程中的边界运动与变形问题，须借助动网格技术来模拟．动网格算法主要包括铺层、弹性光顺及局部重构3种．弹性光顺算法仅限于变形非常小的情况，变形过大时会极度扭曲．局部重构算法适用于大变形或者大位移情况．当网格扭曲过于严重、网格尺寸过小或过大时，网格会自动进行局部重构以满足扭曲和尺寸的要求．针对计算的类型采用弹性光顺和局部重构算法相结合的方法[14]．采用DEFINE_CG_MOTION的UDF函数定义喷针的运动，即指定喷针及折向器的运动速度，并设定喷针的运动为刚体运动．2.3.2　网格变形有效性验证喷针运动属于阀门关闭运动，运动部件与固定边界距离越来越窄，造成间隙越来越小，对应的网格尺寸也越来越小．折向器运动属于普通旋转运动，由于偏转角度较小，对网格引起的弹性变形对计算结果的影响很小，因此须对全局进行网格重构定义，即当网格变形变大到一定尺寸及变小到一定尺寸时，对部分不满足尺寸要求的网格及周边网格进行重新划分，表1为4种网格重构方案．10.13245/j.hust.240653.T001表14种网格重构方案方案最小网格长度/mm最大网格长度/mm最大网格偏度10.325 60.449 4470.933 85320.894 00.449 4470.933 85331.007 00.449 4470.933 85341.325 00.449 4470.933 853对网格重构频率取值为1，保证每一步重构的网格尺寸都能满足计算需求．利用网格变形后，变化的最大网格偏度来辨别网格变形的有效性及质量．网格最大偏度越小，网格质量越高，计算结果越准确；同时网格约束越大，网格重构的次数就越多，每一个步长的计算时间越长．网格约束越大，重构后的网格质量越好，时间越长，根据时间和质量的双重需求，最终选择方案3作为网格重构方案．2.4　一维计算验证法过渡过程计算中采用优化方法进行水轮机水压力迭代计算，先用进退法求出包含解的区域，然后用区域分割法求出解．该方法收敛性好，计算时间短．根据电站输水系统设置的实际情况，建立仿真计算模型．对调压井前的引水隧洞及调压井后的压力钢管进行分段并划分节点，将配水环管作为机组节点接入进行计算．a．有压管道数学模型Qp=Cp-CaHp；Qp=Cn+CaHp；Ca=gM/a ；Cp=Qa+gMaHa-fΔt2DMQaQa；Cn=Qb-gMaHb-fΔt2DMQbQb，式中：Qp和Hp分别为P点当前时刻的流量和水压力；Qa和Ha分别为M点前一时刻的流量和水压力；Qb和Hb分别为B点前一时刻的流量和水压力；A为管道截面积；D为断面水力直径；a为波速；f为摩擦系数．b．有压管道线性模型刚性水击(忽略液体的压缩性与管材的弹性，把管道内的液体视为一条整体的“刚性水柱”，在管道各个截面上的液体流速只是时间的函数，而与位置无关)Gs=Tws+hf ，式中：G(s)为液体的流速；Tw=LVr/(gHr)；Vr为平均流速；Hr为测压管水头；L为管道长度；hf为相对水头损失．弹性水击(考虑液体的压缩性和管材的弹性，在管道各个截面上液体的流速是位置与时间的函数)Gs=b2s2+b1s+b0a3s3+a2s2+a1s+a0；b2=-3/(Trhw)；b1=-6hf/(Tr2hw2)；b0=-24Tr3hw-3hf2Tr3hw3；a3=1；a2=3hf/(Trhw) ；a1=3hf2Tr2hw2+24Tr2；a0=24hfTr3hw+hf3Tr3hw3，式中：Tr=2L/a；hw=aVr/(2gHr)．3 甩负荷三维过渡过程3.1　初始流态分析图4为配水环管初始流态图，从4(a)剖切面流态图可知：配水环管内部及折管内部流速分布较为均匀，配水环管中断面平均流速为13 m/s．10.13245/j.hust.240653.F004图4配水环管初始流态图4(c)为配水环管内壁压力分布图．配水环管内壁所受压力分布基本按照管道延伸方向逐步增大，当流体在通过岔管进入喷嘴时，压力也会逐步递减．根据流量守恒方程及能量守恒方程可解释这一现象，第一岔管末端压力、第二岔管末端压力、第四岔管末端压力均在7.3 MPa以下．3.2　喷针临界关闭工况流态分析图5为喷针开度T1=0.03时各流动参数云图．如图5(a)所示，此时喷针基本完全关闭，入口流量仅为6.184 m3/s，在喷嘴和喷针缝隙出有部分液体泄漏．由图5(b)可知，此时配水环管内部流动基本滞止，管道内流速接近0 m/s．由图5(c)压力分布可知：此时的配水环管管壁所受压力从7.396 MPa至7.493 MPa均匀增加，沿配水环管流向沿程递增，在配水环管末端达到压力最大值，较初始流态的压力值明显降低．此时的剖切面由于流动滞止的原因，其压力分布与管壁压力分布也完全相似．10.13245/j.hust.240653.F005图5配水环管喷针开度T1=0.03时各流动参数云图此开度下，网格最小体积已经到达计算极限，为避免负体积网格出现，此状态即为数值模拟的极限状态．3.3　过渡过程分析配水环管喷针关闭过程为25 s直线关闭，初始开度T1=0.75，至T1=0.1开度截止，因此实际关闭时间为16.25 s，喷针运动速度为14.8 mm/s．监测其运动过程中的流量、压力变化，将压力全部换算为水头(H)来表示．图6为关闭过程中流量(Q)-时间(t)变化曲线图，由图6可知，在关闭过程中流量自78 m3/s均匀减小．折向器关闭0~3s时对流量并无影响，流量变化主要还是由喷针控制．10.13245/j.hust.240653.F006图6关闭过程中流量-时间变化曲线由初期观测可知，管道外壁的压力明显较高，因此将监测点设置在较为危险的部位观测，其位置设置如图7所示．取z=0平面高度的监测点用以表示管道的总体变化趋势．在z=0平面与管道外壁相交的等角度岔管处设置监测点p1~p5，由于配水环管末端无岔管，而是由角度递增的折管起到对水流的导流作用，因此分别在z=0平面与四段折管的外侧相交处取4个监测点p6~p9．对配水环管六出口阵列6个监测点o1~o6观察出口位置压力变化．10.13245/j.hust.240653.F007图7监测点布置图8为管道压力(以水头表示)变化图，由图8(a)~(c)可知：在配水环管各监测点的压力变化趋势基本一致，在折向器动作瞬间均呈下降趋势，当折向器即将完全关闭挡水时，会呈现凸字型的压力分布规律，这一变化在越靠近配水环管末端部位时越发明显．在喷针关闭至0.45开度后，不同部位均会出现持续1~1.5 s的剧烈压力波动，波动还是沿着压力增大的方向变化，到开度0.61附近时，会达到压力上升的第一个峰值，在开度达到0.83附近时，会达到第二个峰值，两处峰值的大小均与监测点的位置及初始压力大小有关系．10.13245/j.hust.240653.F008图8管道压力变化图8(d)为第一岔管处的不同监测点压力变化，其监测点均位于z=0平面．可以发现：岔管外侧压力普遍小于内侧压力，并且产生的波动均较小，在喷针关闭至0.83开度达到第二次峰值后，内外侧的压力逐渐一致，外侧管道压力为一致后内外侧压力差仅为0.5 m．初始压力越低的部位压力上升值越大，上升速度越快．3.4　一维计算结果验证根据电站输水系统设置的实际情况，建立仿真计算模型，将配水环管作为单独节点进行计算．设定计算的初始工况为上游正常蓄水位，两台同时事故甩额定负荷，喷针的初始开度T1=0.754 9，关闭规律为25 s关闭，折向器的初始开度T2=0.985 1，关闭规律为3 s关闭，初始水头为673.66 m，初始流量为83.76 m3/s，与三维计算条件保持一致．将节点的压力绘制成图并与三维计算中各出口处岔管监测点共同排列，结果如图9所示．10.13245/j.hust.240653.F009图9配水环管节点压力由图9可知：整个过渡过程，两种计算其压力遵循相似的变化规律，且对峰值出现的区间基本相对应．相比三维计算结果，一维计算过程更加平滑，对波动的计算精确度较低，而对整体变化的趋势更加敏感．三维计算过程能准确捕捉关闭过程中的压力脉动，但在一维计算中，配水环管被当作一个节点进行计算，其不同部位的瞬态特性很难被具体分析．由此可知，三维计算对于配水环管内部的流态及其特性研究更为合理和准确．4 结论本研究采用欧拉-拉格朗日法对水-气两相流进行了非定常计算．为了验证数模方法的可靠性，对模拟环境的有效性进行了分析，并对网格划分、区域加密控制、网格运动变形尺寸控制、压力等参数变化进行了对比验证，将一维计算加入到对比中．得到的结论如下．a．在水电站大波动过渡配水环管的折向器与喷针同时动作后，其流量依然由喷针控制，折向器的动作不影响配水环管内部的流量变化．b．流速在配水环管内均匀分布，在进入岔管后会因为出现的漩涡导致流速减小，压力则在配水环管内沿管道延伸方向递增，在末端达到最大，并且管道外侧的压力普遍大于管道内侧．在过渡过程后期，配水环管的压力分布趋于均匀，且保持之前的分布规律不变．c．在过渡过程中，在进口段后的所有点压力分布遵循同样的规律，动作开始时出现压力下降趋势，之后会出现凸字型的压力变化规律，在开度关闭至0.45时会出现持续的高频波动，当开度0.61附近会达到压力上升的第一个峰值，在开度达到0.83附近时达到第二个峰值．d．压力变化的幅度与初始压力及部位有关系，在进口段，压力呈逐步上升趋势；在折管段初始压力越大，最大压力越大，且波动越剧烈，压力上升值越小；在岔管段，同一侧的压力会逐渐统一，不被初始压力大小影响．