伊丹河输水河道封冻期冰情演变数值模拟

樊 霖,茅泽育,齐文彪,刘树峰

(1.清华大学水利水电工程系,北京 100084; 2.吉林省水利水电勘测设计研究院,吉林 长春 130021)

伊丹河输水河道封冻期冰情演变数值模拟

樊 霖1,茅泽育1,齐文彪2,刘树峰2

(1.清华大学水利水电工程系,北京 100084; 2.吉林省水利水电勘测设计研究院,吉林 长春 130021)

根据水动力学、热力学、河冰水力学及固体力学等基本理论,针对吉林省伊丹河输水河段具体特征,建立了伊丹河的河冰数值模型,并应用1958—1988共30年的水文、气象资料对该输水河段进行了数值模拟及分析。结果表明,伊丹河输水河段冰期封冻形式为平封,可形成稳定冰盖从而实现冰期冰盖下输水,且满足冰期正常输水的要求。

冰期输水;调水工程;冰情演变;数值模拟; 伊丹河

河冰是寒冷地区尤其是高纬度区域冬季河道中普遍存在的自然现象。北方河道冰期输水时,严寒的气候往往出现冰害,导致运行和管理上的许多困难。因此研究冰期输水所涉及的相关冰害问题对于河渠道冬季安全运行具有重要意义。

河冰生消演变为水动力、机械力和热力等相互耦合作用且极为复杂的物理过程。多年来,随着河冰水力学理论的发展,国内外学者们相继提出了一些模拟河道冰情的数值模型。例如Lal等[1]根据热力和两层输冰模型的理论分析框架,建立了冰情演变河冰数值模型—RICE模型。Shen等[2]提出了二维与一维结合的河冰动力学模型DynaRICE,主要用于模拟冰坝形成的动力学过程等。Zufelt等[3]建立了一维冰水耦合的冰塞动力学模型,模拟封河期和开河期水位变化。茅泽育等[4-5]针对天然河道弯曲复杂的特点,首次建立了适体坐标下的二维河冰数值模型,对黄河河曲段、黄河万家寨水库库区河段、松花江上游白山河段、新疆额尔齐斯河段等进行了数值模拟研究,取得了较好的效果。Fu等[6]在RICE和RICEN模型的基础上,开发了用于模拟冰期洪水的YRIDM模型。靳国厚等[7]采用一维非恒定流水力学及热力学相耦合的模型,预报了京密引水渠、引黄济青等输水工程的冰期输水过程。卞雪军等[8]以多元线性回归理论与人工神经网络理论为基础,设计开发了黄河内蒙古河段冰情预报系统。

如上所述,河流冰情演变受水力、气象、热力、河势及河道特征等众多因素的综合影响,这些影响因素的不同组合必将导致发生不同的冰情。所以,对于任一具体河道,需根据具体的气候、地形、水力等条件建立模型、率定参数并进行数值计算。

伊丹河全流域面积为483 km2,其河源至新立城水库入库断面全长47.71 km(图1)。伊丹河为伊通河的最大支流,沿途流经伊通县二道镇、伊通县伊丹镇、长春市南关区新湖镇等,在长春市南关区新湖镇榆树屯注入伊通河。

图1 伊丹河输水河段示意图

伊丹河冰期输水工程为吉林省中部城市群引松供水长春干线工程的重要组成部分。长春干线全长14.1 km,设计流量13.2 m3/s,由丰满水库经110 km有压隧洞(隧洞取水口位于水库正常蓄水位以下21 m),经冯家岭分水枢纽进入长春干线,再在C09断面处泄入伊丹河,并自流入新立城水库。汇入断面C09至新立城水库入库断面河段长度为7.3 m。由于工程地处严寒地区,河道输水一旦形成冰害,将抬高长春干线末端即出口断面C09的水位,减小长春干线的输水作用水头,影响长春干线设计输水量。由于长春干线与冯家岭分水枢纽连通,还将影响四平泵站、辽源泵站设计参数,并对冯家岭调压井结构产生不利影响,进而影响整个总干线的设计参数。此外,长春干线伊丹河下游末端为水库,极易形成冰害,一旦形成冰害,对河岸、堤防产生影响,影响输水稳定性。因此,实现长春干线伊丹河段冰期安全输水对于长春干线、冯家岭分水枢纽乃至整个中部供水工程正常运行至关重要。本文根据伊丹河河道的地形、气象及研究河段实际情况,建立河冰数值模型,并应用伊丹河1958—1988年的水文、气象资料对伊丹河冰期输水进行了冰情数值模拟。

1 河冰数学模型及计算方法

封冻期河冰演变的模型主要包括河道水力学模型、输冰模型、热力学模型及冰冻模型等。

封冻河道的水流一维流动近似由以下连续性方程和动量方程描述[4]:

(1)

式中:Q为流量;A为过流断面面积;ql为单宽侧向入流量;x为距离;t为时间;g为重力加速度;z为水位;zb为河床高程;ti为冰盖厚度;ρ、ρi分别为水和冰的密度;h为冰盖下水深;nc为综合糙率系数;P为湿周。

根据热量守恒原理,沿河水温η的时空分布近似由以下对流-扩散方程表述[4]:

(3)

其中η=ρCpATw

式中:v为断面平均流速;B0为水面宽;∑S为水体与周围环境的单位面积热交换量(增热或失热);Cp为水的比热容,取4 185 J/(kg·℃);Tw为断面平均水温。

根据分层输冰理论,河流中的流冰由面冰和悬浮冰两部分组成,其质量守恒方程可参见文献[9]。根据冰盖前缘断面的水动力条件,以及上游来冰情况,冰盖的发展一般有并置推进、水力增厚推进、机械增厚推进[1,10]。冰盖下冰的输移和积聚采用文献[11]中的输冰能力公式计算。当积聚冰盖体向上游发展过程中,一方面,冰盖体中的孔隙水冻结而形成固状冰体外壳,增大冰体强度;另一方面,由于冰盖表面和底面发生热交换,冰体厚度将发生热力增厚(或消融)。冰盖体的热力增厚或消融过程由热力平衡原理得到,即:

(4)

式中:ti为水温度；Ts为冰层表面温度;Ta为气温;Tm为冰层热融度,即0℃;hia为大气与冰层表面的热交换系数,取19.71;hwi为水体与冰层底面的热交换系数。

模型中式(1)(2)(3)采用显式有限差分MacCormack步进格式求解[12],其在时间和空间上具有二阶精度o(Δx2,Δt2)。MacCormack方法属两步预测-修正型,其预测和修正步骤,可以交替地向前和向后差分。各子模型组成的联合方程组采用同步求解。为准确模拟伊丹河输水河段的冰情,笔者曾应用白山河段完整详尽的气象、冰情原型观测资料,采用工程类比验证方法对上述数学模型进行了模型参数率定及验证[13]。

自丰满水库取水口至长春干线出水口断面C09之间的埋地输水管道水温,采用以下一维水温对流-扩散方程计算[14],得到长春干线出水口断面水温随时间的变化过程:

(5)

2 伊丹河输水河段数值模拟

2.1 河段基本概况

对伊丹河(全长47.7 km)上游河源、下游新立城水库及丰满水库有压隧洞取水口作全系统研究。长春干线全长14.1 km(PCCP管道,平均埋深6 m),汇入断面C09;丰满水库至长春干线的有压隧洞全长110 km。计算域共布置断面154个，平均断面间距约300 m。断面1布置在伊丹河河源断面(C09断面),编号为84，距河源40.40 km；新立城水库入库断面编号为154，距河源47.7 km；断面127位于河源下游44.93 km。河段平均坡降3.2×10-4。

2.2 气象水文资料选取与整理

采用新立城坝上水文站、伊通气象站1958—1988年水文、气象多年平均观测数据。伊通气象站冬季(11月至次年3月)多年日均气温如图2所示(横坐标时间以11月1日为起点(下同))。计算总时长151 d,时间步长60 s,对于实测断面间距过大河段进行内插,最大计算断面间距取50 m。各实测横断面的过水面积与湿周通过微断面法[15]进行计算。

图2 伊通气象站冬季多年平均日均气温

2.3 初始条件及边界条件确定

初始条件为计算开始时刻计算域内流量、水位、水温、浮冰浓度、水内冰浓度和冰盖体厚度。伊丹河河道初始流量及水位采用上游边界及下游边界对应的流量和水位值,通过长时间迭代运算直至稳定,最终得到沿程流量和水位值;初始水体温度采用6℃。根据水文、气象资料,初始时刻沿程浮冰浓度、水内冰浓度分布及冰盖体厚度均设定为零。流场计算边界条件：上游边界为伊丹河河源断面处的流量过程以及中途汇入的长春干线设计引水流量(13.2 m3/s),下游边界为入库断面(断面154)流量水位关系曲线(图3)。水流温度场计算的边界条件为伊丹河河源断面处以及丰满水库取水口处实测水温值(图4),初始条件通过模型迭代计算,直至初始温度场基本达到热交换的准平衡状态得到。

图3 新立城入库断面冬季流量水位关系线

图4 丰满水库有压隧洞取水口处实测水温随时间变化

2.4 计算结果及分析

图5为封河过程中不同时刻输水河段沿程水位计算结果，图中横坐标以C09断面为起点(下同)。随着气温降低,冰盖前缘逐渐向上游推进,前缘断面附近水位随之壅高,同时冰体厚度增长。1月29日水位壅高达到最大;3月下旬随着气温回暖、水温升高,冰体开始消融,沿程水位降低。图6给出的是C09断面封冻期(11月至次年3月)水位变化过程。由图6可知,12月末1月初,由于下游河段冰盖的生成,C09断面水位开始壅高,最高壅水达0.16 m。随着冰盖进一步固结,冰盖糙率系数降低,计算壅水位逐渐降低。入库断面封冻期流量变化过程如图7所示,12月初,入库断面流量有所减小,这是由于河道槽蓄水量增加所致,随后入库断面流量稳定在13.4 m3/s。入库断面水温随时间的变化过程如图8所示。由于受长春干线输水水体水温的影响,伊丹河输水河段的水温变化滞后于日均气温变化。从图8可知,12月15日受大气温度下降及河道水体与大气热交换的影响,水温持续降至0℃以下,水体进入过冷却状态,随着冰层生成,水温又逐渐稳定在0℃附近。

图5 不同时间沿程水位变化

图6 C09断面冬季水位过程

图7 入库断面冬季流量过程

图8 入库断面水温变化过程

图9 1月29日河道弗劳德数沿程分布

图10 冰盖前缘位置随时间变化

图11 入库断面冰厚随时间变化

图12 1月29日冰厚沿程分布

1月29日输水河段沿程弗劳德数如图9所示。由图9可知,输水河段的弗劳德数均在第二临界弗劳德数(0.09)以下,表明在气象及水温等条件下,河道大部分河段都可形成冰盖。图10给出的是冰盖推进前缘位置随时间的变化，可以看出,12月30日下游入库断面首先形成冰盖,6 d后冰盖前缘到达127断面,由于水温的影响,冰盖前缘到达127断面之后便不在向前发展。图11给出的是封冻期入库断面冰盖平均厚度随时间的变化，可见入库断面的初封时间为12月30日。随后,冰盖体在冰面与大气及冰底面与水体热交换的共同作用下,开始热力增厚。2月20日入库断面冰厚达最大,为0.86 m,此后受气温回暖影响,冰厚开始减小。图12为伊丹河输水河段1月29日冰盖体厚度沿程分布。

3 结 论

a. 伊丹河输水河段冰期沿程弗劳德数均在第二临界弗劳德数(0.09)以下,在气象以及水温条件等适宜条件下,河道大部分河段都可形成冰盖,进而实现河道冰期冰盖下输水,输水河段沿程弗劳德数大多小于0.06,河段封冻形式为平封形式。

b. 河段封冻始于入库断面,然后逐渐上溯,6 d之后冰盖前缘到达127断面,冰期最大冰厚为0.86 m。

c. 在封冻过程中,冰盖前缘推进速度及冰体厚度分布受诸多因素影响,且随冰盖前缘逐步向上游推进,沿程水位也逐渐壅高,最高壅水达0.16 m。冰期入库断面流量值在冰盖形成期间略有下降(河道槽蓄水量增加所致),随后一直维持在13.4 m3/s附近,伊丹河输水河段能够满足冰期设计流量下的正常输水要求。

[1] LAL A M,SHEN H T.A mathematical model for river ice processes[J].Journal of Hydraulic Engineering,1993,117(7): 851-867.

[2] SHEN H T,SU J,LIU L.SPH simulation of river ice dynamics[J].Journal of Computational Physics,2000,165(2): 752-770.

[3] ZUFELT J E,ETTEMA R.Fully coupled model of ice-jam dynamics[J].Journal of Cold Regions Engineering,2000,14(1): 24-41.

[4] 茅泽育,吴剑疆,张磊,等.天然河道冰塞演变发展的数值模拟[J].水科学进展,2003,14(6): 700-705.(MAO Zeyu,WU Jianjiang,ZHANG Lei,et al.Numerical simulation of river ice jam[J].Advance in Water Science,2003,14(6): 700-705.(in Chinese))

[5] 茅泽育,许昕,王爱民,等.基于适体坐标变换的二维河冰模型[J].水科学进展,2008,19(2):214-223.(MAO Zeyu,XU Xin,WANG Aimin,et al.2-D numerical model for river-ice processes based upon body-fitted coordinate[J].Advance in Water Science,2008,19(2):700-705.(in Chinese))

[6] FU C,POPESCU I,WANG C,et al.Challenges in modeling ice floods on the Ningxia-Inner Mongolia reach of the Yellow River,China[J].Hydrology and Earth System Sciences Discussions,2013,10(10): 12293-12329.

[7] 靳国厚,高霈生,吕斌秀.明渠冰情预报的数学模型[J].水利学报,1997,28(10):1-9.(JIN Guohou,GAO Peisheng,LU Binxiu.Ice regime forecast model in open channel[J].Journal of Hydraulic Engineering,1997,28(10):1-9.(in Chinese))

[8] 卞雪军,冀鸿兰,姜新华,等.黄河内蒙古段冰情预报系统的开发[J].水利水电科技进展,2014,34(4):62-65.(BIAN Xuejun,JI Honglan,JIANG Xinhua,et al.Development and application analyze of forecasting system of ice conditions in Inner Mongolia reach of Yellow River[J].Advances in Science and Technology of Water Resources,2014,34(4):62-65.(in Chinese))

[9] 王永填.河冰数值模型及河渠冬季输水冰情特性研究[D].北京:清华大学,1999.

[10] PARISET E,HAUSSER R.Formation and evolution of ice covers on rivers[J].Transaction of Engineering Institute of Canada,1961,5(1): 41-49.

[11] SHEN H T,WANG D S,LAL A M.Numerical simulation of river ice processes[J].Journal of Cold Regions Engineering,1995,9(3): 107-118.

[12] MACCORMACK R W.A numerical method for solving the equations of compressible viscous flow[J].AIAA Journal,1982,20(9): 1275-1281.

[13] 樊霖,茅泽育,吴剑疆，等.松花江白山河段冰情数值计算及分析[J].水科学进展,2016,27(6):890-897.(FAN Lin,MAO Zeyu,WU Jianjiang,et al.Study on the numerical model of river ice for the Baishan reach of Songhua River,China[J].Advance in Water Science,2016,27(6):890-897.(in Chinese))

[14] 崔慧,吴长春.热油管道非稳态工况传热与流动的耦合计算模型[J].中国石油大学学报(自然科学版),2005,29(3):101-105.(CUI Hui,WU Changchun.Heat transfer and flow coupling calculation model of transient scenario for hot oil pipeline[J].Journal of China University of Petroleum,2005,29(3):101-105.(in Chinese))

[15] FAN L,MAO Z Y,BAO J,et al.A river-ice model for Wanjiazhai reach of the Yellow River[C]//Proceedings of the 22nd IAHR International Symposium on Ice.Singapore:NUS,2014:869-876.

Numericalsimulationofriver-iceprocessesinawater-transferriverofYidanRiverduringfreezingperiod//

FAN Lin1, MAO Zeyu1, QI Wenbiao2, LIU Shufeng2

(1.DepartmentofHydraulicEngineering,TsinghuaUniversity,Beijing100084,China; 2.JilinInvestigationandDesignInstituteofWaterResourcesandHydropower,Changchun130021,China)

Based upon the theory of hydrodynamics, river mechanics, thermodynamics, river-ice hydraulics and mechanics of solids etc., a mathematic model of the river-ice processes for the Yidan reach is constructed. 30-year long-term series comprehensive field observation data, including geomorphy, meteorology, hydrology, thermo-dynamics and ice regime, are used to simulate the water transfer process of the Yidan River during freezing period. The results show that steady ice cover can be formed by juxtaposition on the Yidan River, which achieves the goal of water transfer beneath the ice cover. In addition, the design discharge of the Yidan River water transfer project can be guaranteed during freezing period.

water transfer during freezing period; water diversion project;river ice process; numerical simulation; Yidan River

国家重点研发计划(2016YFC0402504);清华大学博士生短期出国访学基金(2013118)

樊霖(1990—),男,博士研究生,主要从事河冰水力学研究。E-mail:l-fan11@mails.tsinghua.edu.cn

茅泽育(1962—),男,教授,博士,主要从事水力学及河流动力学研究。E-mail:maozeyu@tsinghua.edu.cn

10.3880/j.issn.1006-7647.2017.06.003

TV143

A

1006-7647(2017)06-0014-05

2017-02-21 编辑：郑孝宇)