冰期明渠非恒定流数值模型
1闫弈博，1冷星火 黄会勇，2王长德
（1长江勘测规划设计研究院，武汉 430010；2武汉大学水利水电学院，武汉430072）
作者简介：闫弈博 （1986-），河南洛阳人，工程师，博士，主要从事水力学文学的研究。
摘要：在输水渠道明渠过水断面上的流量、流速及水位等水力要素随时间不断变化，形成明渠非恒定流流态，而以质量守恒及动量守恒定律为基础的Saint-Venant，是百年来描述一维明渠非恒定流主要手段。对于冰期渠道而言，随着冰盖的形成和发展，湿周及水力半径同敞流期相比均发生变化；本文以圣维南方程为基础，研究了冰盖发展过程中明渠非恒定流数值模型。
关键词：冰期 明渠 水力学 非恒定流 数值模型
1 引言
南水北调中线工程总干渠从北纬32°的丹江口水库引水，跨越8个纬度，穿越长江、黄河、淮河、海河四大流域，途经豫、冀、京、津四个省（直辖市），将水输往位于北纬40°的北京市，属特大型跨流域调水工程。中线输水线路自南向北由暖温带流入半寒冷地区，冬季运行时受寒冷气温条件影响，将有不同程度的冰情出现，对渠道运行安全、运行效率造成影响，对渠道控制运行提出了新的要求。
冬季随着气温逐渐降低，冰盖形成后造成明渠可过水断面减小，挤占了渠池可容纳蓄量；另外，随着浮冰盖形成，水流湿周增大、水力半径减小；故即使没有冰坝产生，渠道水位也会有所抬高。若冰盖下运行期渠道中发生流量、水位变化，可能对冰盖安全造成威胁，故需要针对冰期水力条件整理修正圣·维南方程。
2 水流阻力分析
河流力学对水流阻力的划分包括：河床阻力、边壁阻力、河床形变阻力（弯道、收缩等）、床沙运动阻力。对于冰期河渠，还包括风、冰盖、岸冰、锚冰等造成的阻力。
对于明渠输水而言，一般为混凝土衬砌渠道，在水源为水库水的情况下，可以不考虑床沙运动阻力；另外由于输水明渠中多有倒虹吸、节制闸等建筑物，因此明渠水流阻力除了床面阻力、边壁阻力及河床形变阻力三项外，还应包括建筑物阻力。
河渠非恒定流运动方程一般为：
当渠道内冰凌含量较少时冰期渠道水力条件受冰凌影响很小［2］，主要考虑冰盖下输水过程水力条件的变化。
忽略冰的相态变化对水体质量守恒的影响，对于具有浮动冰盖的河流，其水位和流量的变化可以由连续性方程和动量方程描述。
式中：τa为水体表面的风拖曳力；τb为水流对河床的剪切应力；τi为水流对冰的剪切应力；pi为冰盖湿周；pb为河床湿周；ρw为水体密度；B0为不包括岸冰的水面宽度。
很多学者都研究了风对冰面的拖曳力τa，但是相对来说，风拖曳力的量级一般为10-2N/m2，可将其忽略［2］。
3 冰封段水流阻力
3.1 阻力计算公式
式（1）中，水力坡度J是由水流受到的各种阻力造成的。如果可以通过实验或者现场测定得到各个阻力项，那么运动方程就可以成功建立并进行求解。但是，由于阻力影响因素的时空变异性及不确定性，这常常是不可能完成的任务。为了建立可以求解的运动方程，引入了基本假设［1］：
（1）假定过水断面符合静水压强分布规律。
（2）假定局部水头损失可以忽略，沿程损失可以用恒定流阻力公式计算。
（3）假定过水断面上流速均匀分布。
（4）假定河床为定床，且底坡≤0.1。
在以上四条假设基础上，可以采用恒定流阻力公式计算非恒定流水力坡度。
国内外学者建立了多种公式计算恒定流水力坡度，常用的包括达西-魏斯巴哈公式及谢才公式。
（1）达西-魏斯巴哈公式：
（2）谢才公式：
谢才公式在我国明渠水力学中得到了广泛的应用，因此本文也采用谢才公式计算水力坡度，并采用曼宁公式计算谢才系数：
曼宁公式：
将式（5）、式（6）联立整理，可以得到：
3.2 综合糙率
在式（7）中，糙率n代表了水流受到的全部阻力，故又可以称为综合糙率系数。
孙肇初、隋觉义等［3］根据黄河河曲段冰期实测资料，得到河段冰期综合糙率的经验公式：
Sabaneev在将流动分区分为床面影响区、和冰盖影响区的基础上，假定过流断面各处流速相等［4］，得到综合糙率公式：
式中：ni、nb分别为冰盖下表面及河床糙率，其他变量意义同式（3）。
在对多个综合糙率公式进行分析比较的基础上，考虑不完全封冻情况，如图1所示，河流中部分封冻河流断面可以划分为两个部分。
图1 河冰横剖面图
魏良琰［5］引入封冻系数fb=B0/B，并引入基本假定。
（1）过流断面随湿周糙率的不同可分为几个相应的分区，以最大流速线为分界。
（2）因为二次流的作用，相同糙率的表面上具有相同的均布剪应力。
（3）均匀流阻力规律同样适用于恒定、非恒定流。
（4）分区水力坡度相等并等于总流水力坡度。
（5）分区平均流速相等并等于总流平均流速。
在以上假定的基础上，改进了Sabaneev公式：
对于输水明渠而言，由于断面尺寸、沿程底坡变化均较小，弯道较缓。因此，假定冰盖仅在拦冰索等拦冰建筑物前发端，并且冰盖发展速度没有横向差异，即有冰封渠段fb=1。则对于冰盖下明渠非恒定流而言，综合糙率计算的关键是确定渠冰覆盖范围。
4 冰盖发展段水流阻力
4.1 明渠冰盖范围
对于长距离输水明渠而言，干渠可认为被节制闸及建筑物前的拦冰索划分为若干个冰盖发展区段，这些区段的冰盖发展相互独立；而冰盖下非恒定流模拟以可测量、控制的渠池内上下游边界划分计算段，因此圣维南方程计算段内可能同时存在敞流段及冰封段，故其引用综合糙率，还需要进一步考虑渠段内冰封范围对综合糙率的影响。
在一个区段中，冰盖在拦冰索等拦冰建筑物前发端，整个区段被划分为两部分，敞流段及冰封段，如图2所示。
在进行水力学仿真时，对于一组边界条件之间，长度为ΔX的渠段：设圣维南方程阻力系数、或者说水力坡度为一恒定值J；而冰盖推进是以特征线距离步长为单位；因此ΔX长的渠段内有若干距离步长为冰封段，其他为敞流段。
图2 封河过程
图2中，区段长度为ΔX，若有敞流段水力坡度Jw，及冰封段水力坡度Ji，敞流段占渠段比例为λ，则区段平均水力坡度可以求得：
在明确冰盖覆盖范围的影响后，还需要确定冰盖糙率。
4.2 冰盖糙率
冰盖在初步形成后糙率较大，之后随着水流的摩擦而变小并趋于一个常值。在开河前会因融冰而造成冰盖下表面糙率增大。冰盖糙率需要通过原型观测得到，Nezhikhovskiy［6］在分析研究苏联大量原型观测资料的基础上，提出了冰盖糙率的计算方法：
初封期冰盖，ni=0.01～0.02。
开河前冰盖，ni=0.008～0.01。
在此基础上，假定稳定冰盖期的冰盖糙率变化为一个指数衰减过程：
式中：ni，e为冬季末期的稳定冰盖糙率；ni，i为冰盖形成初期的冰盖糙率；k为衰减系数；t为封冻天数。
式（12）的缺点在于ni，i和ni，e难以确定，并且单一的指数k难以描述整个冰期冰情的剧烈变化。
孙肇初、隋觉义［3］根据黄河河曲段原型观测资料，得到经验公式：
式（13）建立了冰盖糙率同河床糙率、断面平均流速及水深的关系，但是参数的率定随着河床变化将会是一个难题。
沈洪道等［7］根据St.Lawrence河的原型观测资料，分析认为冰盖阻力决定于水流条件、冰盖情况及气候条件等因素。在冰盖糙率先升后降的整个变化过程中，会因气候、冰花输移的变化而产生脉动，即
式中：nd为基础冰盖阻力，封冻期单调增加、稳封期单调减小；nt为融化期的冰盖阻力增量；n′为冬季水内冰输移和沉积的一个波动分量。
5 冰封渠道圣维南方程
结合上述分析，将式（11）带入冰封河流阻力项，可以表示为
其中冰封河段水力坡度、敞流段水力坡度都由式（7）计算。
对于冰封段糙率：
敞流段糙率：
则，式（12）、式（16）、式（17）及式（7）联立计算并带入式（15），即得到统一格式的冰封渠段圣维南方程。
6 结论
本文从渠道冰期运行水力学模拟的需要出发，分析研究了冰封渠段及冰盖发展段水流阻力，并进一步推导了能够反映冰盖发展过程的圣维南方程统一形式，为冰期渠道水力学模拟提供了模型基础。
参考文献
［1］ 魏良琰.封冻河流阻力研究现状［J］.武汉大学学报，2002（2）：1-9.
［2］ Shen Hungtao,Wang Desheng.Under Cover Transport and Accumulation of Frazil Granules［J］.Journal of Hydraulic Engineering,1995.121(2):184-194.
［3］ 孙肇初，隋觉义，等.水内冰冰塞形成和演变的实验研究［R］.国家自然科学基金资助项目研究报告，1990.
［4］ WangJun.Discussion of velocity distribution and resistance under ice cover.Proceedings of the 6th international Conference on hydrology.Hebei,China.Jun2004:367-370.
［5］ 魏良琰，黄继忠.冰盖流阻力与综合Manning糙率［J］.武汉大学学报，2002（8）：1-8.
［6］ Nezhikhovskiy，R.A.1964.Coefficients of roughness of bottom surface on slush-ice cover.Soviet Hydrology,Washington,Am.Geoph.Union,127-150.
［7］ Shen H T,Yapa P D.Flow Resistance of a River Ice Cover［J］.Journal of Hydraulic Engineering,ASCE,1985,112(2):142-156.