基于离散元方法的河冰动力过程短期预测模型及应用

翟必垚 张宝森 季顺迎 田治宗

摘 要：秋末冬初时黄河河道平稳封河对于预防凌汛意义重大，而什四份子作为黄河内蒙古河段主要首封河段备受重视。根据寒区河道具体冰情，基于河冰离散元方法，考虑河流表面冰的动力输移和河流水动力学的相互耦合作用，建立了针对特定重点河段的短期预测模型，可根据实际水力要素、气温、现场冰情等条件进行河冰输移、封河开河等情况的模拟和预测。模型考虑了河道的岸冰边界以及上游来冰情况，包括来冰的密集度、冰块体的尺寸分布、冰块体的冻结强度等因素对河冰输移和卡冰封河过程的影响。基于该预测模型，针对什四份子弯道在不同工况下的封河情况进行研究，结合现场情况探究了流速、来冰量、冰块尺寸、岸冰尺寸等因素对卡冰封河的影响，从冰流量和河道输冰量的角度解释了各要素的影响机理，为特定河段的封河预测工作提供参考。

关键词：河冰;封河;离散元方法;数值模型;什四份子弯道

中图分类号：TV133;TV882.1

文献标志码：A

doi：10.3969/j.issn.1000-1379.2021.10.015

引用格式：翟必垚，张宝森，季顺迎，等.基于离散元方法的河冰动力过程短期预测模型及应用[J].人民黄河，2021，43（10）：75-80.

Abstract： Based on the discrete element method， a short-term prediction model of the nature river ice dynamic process was established， considering the coupling effect between surface ice and the river hydrodynamics. It could simulate and predict the several river ice dynamic processes， such as， river ice transport， ice bridge formation， ice jam， break up and jam release. Based on the prediction model， the formation of ice bridge in Shensifenzi reach was simulated and the effects of flow velocity， ice supplement， ice size and border ice size on the jamming were investigated. The mechanism of each factor was explained by ice discharge and ice transport capacity of the channel. This paper provided a reference for the forecast work of ice bridge formation and river freeze up in natural rivers.

Key words： river ice; freeze up; discrete element method; numerical model; Shensifenzi reach

冰凌在我國北方寒冷地区的江河非常常见，是气象水文、河道地形地势和水力条件综合影响的结果。严重的冰情会直接影响江河的航运交通、灌溉供水、水力发电，带来巨大的经济损失，还会诱发凌汛灾害，对人民生命财产造成巨大的损失[1-3]。秋末冬初时节，气温降低，江河出现流凌，随着低温的持续，会发生封河现象。一旦封河，应及时减小流量争取平封，防止立封和冰塞的产生[4]。什四份子弯道是黄河内蒙古河段高频的首封河段，其地势形态特征典型，上宽下窄，弯变急且角度变化大，是防凌重点关注的位置[5-6]。该河段已安装配套了多种监测设备进行冰情、气候、水力要素等的实时监测，并配合无人机、卫星雷达定期观察冰情发展，为冰情预测提供了充足的参考依据[7-9]。

国内外的学者基于冰力学理论基础开发了多种河冰热力学、动力学过程模拟的数学模型和数值模型，较为系统的有RICE[10]、DynRICE[11-13]、RIVJAM[14]、ICEJAM[15] 、ICESIM[16]、RIVICE[17]、River1D[18]、HEC-RAS[19]等模型。已有的河冰模型从河冰动力学、热力学机理出发，能够很好地模拟长时段、长距离的河冰热力生消、动力输移等演化过程，对冰塞冰坝的出现位置及引起的水位、流量变化有较好的模拟结果[20-22]。近年来随着人工智能技术的发展，统计模型、人工神经网络、模糊数学等方法也在冰情预测预报中得到了应用，并取得了很好的效果[23-25]。然而，以上模型和方法在实际应用中可能会受到长距离、多时段监测数据不完备的限制，且关注点多集中在大中尺度的水位、流量和冰塞位置上，针对特定局部河段的冰情精细化模拟和预测相对较少。

黄河防凌监测技术发展迅速，已由早年的监测手段少、监测指标单一发展形成了初步的冰情观测技术体系，测量指标包含流速、水深、冰厚、水面冰密集度、水温、冰温等[21]，监测手段扩展到了卫星、雷达、无人机、远程4G视频拍摄等[26]。这些都为精细化预测模型的参数设置提供了依据。

本研究基于河冰离散元模型，建立面向真实河道的数值模型，根据特定天然河道具体冰情，对封河前期的水力要素、地势地形、初始冰情等重要条件建立模拟参数，形成针对小范围特定重点河段的短期预测模型;结合什四份子险工处已建立的完备监测手段，对该河段河冰输移及封河情况进行模拟研究，以期为防凌预报提供有效的参考。

1 基于河冰离散元方法的河冰动力过程短期

预测模型

1.1 河冰动力学短期预测模型的计算框架

结合特定河段的监测数据，基于河冰离散元方法，建立特定河道的短期预测模型，可对特定局部河段的河冰动力学过程进行精细化模拟，具体模型框架如图1所示。

由于热力计算的输入条件较为复杂，涉及日照、气温、水温等条件，且热力计算多针对长时间尺度的模拟、误差不易控制，因此模型回避了热力学计算，将河冰的热力学生效结果转化为动力学模拟的初始条件，如岸冰边界、冰块冻结强度等，再通过短期的动力学模拟以提高模型的适用性和准确性。完备的河段冰情资料包含气象、水文、河道地形、冰情等條件，其可对模型的正确性进行验证，也可对模型计算的初始参数进行校核，进而对特定的河冰现象如封河、开河、冰坝、溃坝等进行模拟计算，为河冰现象的预报提供参考。

1.2 河冰动力过程的离散元方法

河冰动力过程计算模型考虑了河冰动力学与水动力学的耦合作用[27]。河冰动力学采用了二维河冰输移模型，可计算河冰受到的水流和风的拖曳力、水面梯度力、冰块内力、冰与河岸及建筑物的接触力;河流水力学采用Shen根据表面冰层作用改进的二维浅水模型[11];河冰内力主要采用离散元方法显式逐步求解[28]。此外，该模型考虑了河冰离散单元的浮力、拖曳力以及水流与冰层的耦合[27]。

根据河冰在低温条件下的冻结现象，建立了河冰离散元黏聚力模型，即计算单元之间建立黏结对，使得单元和单元之间可以相互黏结在一起，黏结对可以在单元之间传递拉力、压力、剪切力、弯矩和扭矩，当传递的力超过破坏准则时，黏结对就会失效，宏观表现为冻结冰的破碎。扩展多面体单元的黏聚力模型有两种类型（见图2）。一种是建立在共面接触的黏结单元间，即必须是面面接触模式才可建立黏结对[28-29]，黏结面积即为接触面积。另一种是针对自然状态下冰堆积体中的冰块接触随机不规则的情况，建立的基于任意接触方式离散元黏聚力模型，这种黏聚力模型在单元之间建立虚拟黏结面，虚拟黏结面的面积直接影响黏结对传递的力的大小，其主要受到气温T和累计接触时间ta控制，其表征了冰块由接触到冻结以及冻结面生长的过程。具体关系需要根据特定河道的现场试验结果或者相关规范进行选取。

根据河冰在动力过程中表现出的离散分布特性，采用Voronoi切割算法结合扩展多面体离散单元形态对冰块体进行构造。二维Voronoi切割算法是在一个二维平面内随机撒点，再通过随机点将该二维平面划分为随机的任意多面形，对多面形赋予冰厚便可得到随机分布的碎冰模型。碎冰的尺寸、形状以及碎冰分布的密集度可以根据随机点的分布参数进行调节[30]。

1.3 河冰动力过程离散元模拟的边界条件

初始岸冰的构造主要依据现场观测的岸冰形态，通过绘制岸冰边界，再采用随机三角单元对岸冰区域进行剖分，组成岸冰的三角冰块初始状态是相互黏结的，由此得到岸冰模型。在岸冰与河岸接触的位置采用固接的方式将二者黏结。由于岸冰的热力生成过程长、冻结时间久，通常有较大的强度，因此岸冰区域的冰单元黏结强度通常设置为较大值。

根据不同算例工况对冰块补给方式的需求，模型中给出了两种冰块补给方法，分别为周期边界法和上游供给法[31]。周期边界法能够大大节省程序计算内存以及消除边界对颗粒流动的影响，被广泛应用于离散元模拟中。周期边界法可应用在规则河道的两岸以消除边壁效应，也可应用于河道上下游边界处，以节省程序预构造冰块数量。这里主要针对后者的应用进行说明，其实质为当冰单元从河道下游边界流出后，会从上游重新进入河道区域。图3为周期边界法中的河冰输移。为了保证边界处的冰单元可以正确地搜索并接触到另一端边界处的冰单元，根据边界附近的冰单元（Z1，Z2，Z3，Z4）在另一边界处建立相应的虚拟冰单元（V1，V2，V3，V4）。如果河道为规则直河道，则可通过平移实现;当河道上下游边界方向不一致时，可采用四元数控制旋转角度的方法通过平移和旋转来实现。

上游供给法是另一种用来补给河道冰量的方法，即在河道的上游边界或者特定位置进行输冰。程序初始化阶段预先构造足够多的冰块单元，开辟相应内存但不参与计算。在河道上游边界或者特定位置设置来冰边界（见图4），确定来冰边界输冰个数Nb和输冰间隔dl0。边界输冰个数Nb是指输冰边界处一个时间间隔最多能够输送Nb块碎冰，冰块的输冰间距dB0 =B/Nb;输冰间隔dl0可以控制当输送的冰块运动多远距离后开始进行下一轮输冰。由于在输入冰单元时必须保证输入冰块之间以及输入冰块与河道已有冰块之间没有接触重叠，这里采用了最大颗粒球形包围盒尺寸Dmax控制冰块间隔，即dB0

周期边界法多用于规则河道的算例研究，而上游供给法冰块由上游边界输送，更符合实际情况，在自然河道算例中推荐采用。不论采用哪种冰块补给方式，计算中都会根据初始时刻实际河道的冰分布情况，预先在河道内布置一定量的冰块，这种方式可减少冰块由上游边界流动到占据整个河道的时间，大大提高了计算效率。

2 黄河内蒙古什四份子河段封河过程的离散元模拟

2.1 黄河内蒙古什四份子河段概况

什四份子河段地处黄河内蒙古河段下游（北纬40°17′39″，东经111°2′53″），见图5。该河段地形特征明显，由上游的宽河道在入弯处逐渐束窄，河道走势方向也发生了剧烈变化，因此该位置较易产生卡冰现象，历史上多次成为秋冬封河期首封位置。

根据2019—2020年该河段封河过程的观测资料，2019年什四份子封河时间为12月6日，封河流量约为787 m3/s，平均流速约为0.71 m/s，封河当日气温为1.4～-15.3 ℃，平均气温为-10.1 ℃;而2020年封河时间为12月8日，封河流量约为362 m3/s，当日气温为-7.2～-18.2 ℃，平均气温为-13.6 ℃，平均流速约为0.38 m/s。这两年封河当日冰况差别较大，图6为什四份子河段2019年和2020年封河过程的现场航拍图。2019年什四份子河段封河之前上游来冰量较大、冰块尺寸小且密集，河岸边已有一定数量的岸冰生成;2020年封河之前来冰呈大块片状，冰块稀疏但是尺寸较大，在持续低温的作用下岸冰生长也更加明显，与2019年情况相比更加束窄了河道。

2.2 短期封河的离散元数值预测

根据2019年和2020年封河前一天的水力情况和冰情，建立模型模拟封河现象，模拟参数见表1。模型中水力边界条件上游为流量，下游为水位，平均糙率0.03。图7为计算网格及河床高程。采用恒定流研究河冰的输移特性，2019年工况设定初始流量787 m3/s，下游水位988.2 m，稳定流态如图8（a）所示，主槽的平均流速为0.7 m/s，上游输入冰块的平均面积为500 m2，厚度为0.2 m，密集度为0.8左右;2020年工况将初始流量调整为362 m3/s，得到的稳定流态如图8（b）所示，主槽平均流速约为0.4 m/s，输入冰块平均面积为5 000 m2，厚度为0.2 m，密集度为0.3左右。

图9为这两年封河模拟的结果。可见在相应水力条件和来冰量的情况下，冰块在河道束窄位置卡冰，模拟结果与实际观测现象（见图6）较为吻合。从结果中可以看出岸冰的束窄作用决定了卡冰的位置并直接导致卡冰的发生，冰块大小和流凌密度也影响了卡冰难易程度。2019年什四份子河段是由小尺寸高密度来冰在大流量下发生卡冰封河的，而2020年该河段是在小流量低流速的情况下因大尺寸的来冰以及岸冰对河道的束窄作用最终导致卡冰封河的。

图10为2019年封河模拟结果中统计的弯道出口处（图9中红虚线所示）冰流量变化情况。在流凌初期河道将来冰顺畅地输往下游，随着来冰量的持续增加，冰流量随之增加并到达稳定值附近，该稳定值即反映了河道的输冰能力。在1.0～2.3 h内，主流继续保持最大输冰能力输冰，出弯处冰流量基本不变，而超过河道输冰能力的冰块会在主流两侧停留、聚集，最终导致河冰在束窄位置出现拥堵。一旦出现卡堵的情况，河道输冰能力会迅速降低，冰流量随之减小，进而卡堵封河。由图10可知不足0.5 h时间内出弯处的冰流量由稳定值降为0。2020年的冰流量变化情况与2019年较为类似，这里不做赘述。

2.3 河冰封河过程的主要影响因素分析

基于对实际工况的分析得知，影响封河的因素较为复杂，包括水力要素（主要是水流流速）、冰量、冰尺寸以及岸冰的边界约束等。为了理解卡冰封河潜在的物理机制，通过数值模型对更多工况进行模拟，进而对封河的条件建立一定认识。

首先考虑无岸冰的工况。图11是密集度为0.4和0.8的来冰输移情况，分别代表低密度来冰量和高密度来冰量。当密集度较小时，冰块主要在河道主槽输移，只有极少量的碎冰停留在岸边低速区。当密集度增大时，主槽的输冰能力无法输送大量来冰，会有部分冰块被排挤至主槽两侧的低流速区域，在弯道凹岸处逗留，并且随着来冰的持续增加，凹岸的冰块不断聚集，无法及时流出。二者都没有出现卡冰封河现象。

当主槽的平均流速降到0.4 m/s左右时，输移情况没有太大变化，河道仍没有发生卡冰封河现象。由于流速降低，主槽的輸冰能力随之降低，因此有更多的冰块在河道两侧停留。

在无岸冰的几种工况下都未发生卡冰封河，可见岸冰的生长对卡冰封河有必要性。流凌初期在持续输冰下，会有部分冰块在凹岸侧停留，来冰量越大，进入凹岸低速区的冰块越多，这给什四份子弯道凹岸处冲击岸冰和平滑岸冰的生长提供了基础，进而为封河的发生提供了有利条件。由此得出：在流凌初期该处河道卡冰封河的可能性较小，这时应主要关注岸冰的生长发展趋势。

其次，根据2019年、2020年两种岸冰形态，探究在岸冰存在的条件下，不同流速和来冰量的封河情况。结果显示，在同样的条件下，密集度越大、流速越慢时，越容易卡冰。这里统计了两种岸冰状态下0.1、0.4、0.8、1.2 m/s几种流速封河时的临界来冰密集度，见图12。较高的流速需要较高的临界来冰密集度才能满足卡冰封河的条件，并且岸冰二较岸冰一的尺寸更大，其对河道的束窄作用更明显，因而对应的封河临界来冰密集度较小，更易形成卡冰封河。

3 结 论

针对天然河道及工程条件，建立了基于河冰离散元方法的短期河冰预测模型，可准确表征实际河道的地形地势、水力要素、冰况信息。将该方法应用于黄河内蒙古段什四份子弯道处河冰输移的模拟和研究中，结合2019年、2020年两年封河前的工况条件进行封河预报模拟，得到的封河结果和卡冰位置与实际情况吻合较好。进一步对该河段的封河条件进行研究，建立了两种特定岸冰形态和特定冰尺寸下不同水流速度对应封河临界来冰密集度的关系，得到了一些定性的结论：河道的边界直接影响卡冰封河的可能性，什四份子河段为由上游宽大到下游狭窄的形态，流冰容易在弯道凹岸的低速区停留，形成岸冰。岸冰的形成是卡冰的必要条件之一，岸冰对河道束窄得越严重，越容易造成卡冰;大尺寸的冰块加速了卡冰封河的过程;流速决定河道的输冰能力并影响卡冰的稳定性，高流速下河道输冰能力强，冰块会沿主槽向下游输移，或直接被带入水底、流往下游，不易稳定封河，低流速下河道输冰能力弱，冰块容易在束窄处聚集形成堆积之势，进而稳定，拦截来冰形成封河。本研究为河冰现象短期预测提供了新的方法，对封河条件的研究成果也可以为封河预警工作提供一定的参考。

参考文献：

[1] 杨开林.河渠冰水力学、冰情观测与预报研究进展[J].水利学报，2018，49（1）：81-91.

[2] 王军，陈胖胖，江涛，等.冰盖下冰塞堆积的数值模拟[J].水利学报，2009，40（3）：348-354.

[3] BELTAOS S. Hydrodynamics of Storage Release During River Ice Breakup[J]. Cold Regions Science & Technology， 2017，139： 36-50.

[4] 蔡琳.中国江河冰凌[M].郑州：黄河水利出版社，2008：5-7.

[5] 罗红春，冀鸿兰，郜国明，等.黄河什四份子弯道冰期水流及冰塞特征研究[J].水利学报，2020，51（9）：1089-1100.

[6] 趙水霞，李畅游，李超，等.黄河什四份子弯道河冰生消及冰塞形成过程分析[J].水利学报，2017，48（3）：351-358.

[7] 武桂芝，张宝森， 李春江，等.阵列地质雷达在黄河堤防隐患探测中的应用[J].人民黄河，2020，42（8）：113-116.

[8] 刘辉，冀鸿兰，牟献友，等.基于无人机载雷达技术的黄河冰厚监测试验[J].南水北调与水利科技（中英文），2020，18（3）：217-224.

[9] 张磊，牟献友，冀鸿兰，等.基于多波段遥感数据的库区水深反演研究[J].水利学报，2018，49（5）：639-647.

[10] LAL W，SHEN H T. A Mathematical Model for River Ice Processes[J].Journal of Hydraulic Engineering， 1991，117（7）： 851-867.

[11] SHEN H T， SU J， LIU L. SPH Simulation of River Ice Dynamics[J].Journal of Computational Physics， 2000，165（2）： 752-770.

[12] SHEN H T ， WANG D S， LAL W. Numerical Simulation of River Ice Processes[J].Journal of Cold Regions Engineering， 1995， 9（3）：107-118.

[13] WANG D S， SHEN H T， CRISSMAN R D. Simulation and Analysis of Upper Niagara River Ice Jam Conditions[J].Journal of Cold Regions Engineering， 1995， 9（3）：119-134.

[14] BELTAOS S. Numerical Computation of River Ice Jams[J]. Canadian Journal of Civil Engineering， 1993，20（1）： 88-99.

[15] FLATO G M， GERARD R. Calculation of Ice Jam Thickness Profiles[C]//Proceedings of 4th Workshop on Hydraulics of River Ice. Montreal， Quebec， Canada， 1986（C3）：1-25.

[16] CARSON R， BELTAOS S， GROENEVELD J， et al. Comparative Testing of Numerical Models of River Ice Jams[J].Canadian Journal of Civil Engineering， 2011，38（6）：669-678.

[17] LINDENSCHMIDT K， XU J Y. Rivice-A Non-Proprietary， Open-Source， One-Dimensional River-Ice Model[J].Water， 2017， 9（5）：314.

[18] BLACKBURN J， SHE Y. A Comprehensive Public-Domain River Ice Process Model and Its Application to a Complex Natural River[J].Cold Regions Science and Technology， 2019， 163：44-58.

[19] BRUNNER G W. HEC-RAS（River Analysis System）[C]//North American Water and Environment Congress & Destructive Water. Washington D C： ASCE， 2002：3782-3787.

[20] SHEN H T. Mathematical Modeling of River Ice Processes[J].Cold Regions Science and Technology， 2010. 62（1）： 3-13.

[21] 郭新蕾，王涛，付辉.河渠冰水力学研究进展和动向[J].力学学报，2021，53（3）：655-671.

[22] BELTAOS S. Progress in the Study and Management of River Ice Jams[J].Cold Regions Science and Technology， 2008， 51（1）：2-19.

[23] GUO X， WANG T， FU H， et al. Ice-Jam Forecasting During River Breakup Based on Neural Network Theory[J].Journal of Cold Regions Engineering， 2018， 32（3）：1-8.

[24] WANG T，YANG K，GUO Y.Application of Artificial Neural Networks to Forecasting Ice Conditions of the Yellow River in the Inner Mongolia Reach[J].Journal of Hydrologic Engineering， 2008， 13（9）：811-816.

[25] MORSE B， HESSAMI M，BOUREL C. Mapping Environmental Conditions in the St. Lawrence River onto Ice Parameters Using Artificial Neural Networks to Predict Ice Jams[J].Canadian Journal of Civil Engineering， 2003， 30（4）：758-765.

[26] 张宝森，张兴红.冰凌监测新技术研究进展[J].中国防汛抗旱，2017，27（6）：14-19，33.

[27] 翟必垚，刘璐，张宝森，等.基于离散元方法与水动力学耦合的河冰动力学模型[J].水利学报，2020，51（5）：617-630.

[28] 刘璐，尹振宇，季顺迎.船舶与海洋平台结构冰载荷的高性能扩展多面体离散元方法[J].力学学报，2019，51（6）：1720-1739.

[29] LIU L， JI S. Bond and Fracture Model in Dilated Polyhedral DEM and Its Application to Simulate Breakage of Brittle Materials[J].Granular Matter， 2019，21（3）：1-16.

[30] SOTOMAYOR O E，TIPPUR H V. Role of Cell Regularity and Relative Density on Elasto-Plastic Compression Response of Random Honeycombs Generated Using Voronoi Diagrams [J].International Journal of Solids and Structures， 2014，51（21-22）：3776-3786.

[31] BABIC M. Discrete Element Simulations of River Ice Transport[C]//Proceedings of IAHR the 10th International Symposium on Ice. Espoo， Finland： IAHR，1990：564-574.

【責任编辑 许立新】