河流水冰沙耦合模型研究Ⅱ：验证与应用

潘佳佳，Hung Tao Shen，郭新蕾

（1.流域水循环模拟与调控国家重点实验室，中国水利水电科学研究院，北京 100038；2.Department of Civil and Environmental Engineering，Clarkson University，Potsdam，NY 13699-5710）

1 研究背景

北方河流冬季流凌能引起严重岸滩刮擦侵蚀和崩塌破坏，伴随着河冰运动的洪水威胁着河道堤防安全。例如，我国黄河和黑龙江等河流的凌汛时有发生，给人民生命和财产安全带来巨大隐患［1-3］。黄河部分河段遇到不利气温和大水条件时，流凌或冰坝能造成堤防溃决，进而引起漫堤洪水，形成严重凌灾［4-7］。此外，流凌刮擦不仅影响岸滩稳定性，也能破坏水工建筑物，危及沿岸建筑或房屋安全。研究冬季河冰运动和岸滩刮擦侵蚀的作用机制对防凌减灾具有重要意义。

国内外很少有研究关注河冰运动对岸滩的刮擦侵蚀机制，已有的部分研究仅定性分析了河冰运动规律或岸滩侵蚀现象［8-9］。Lawler 在南威尔士伊尔斯顿河的两个弯道上开展了长达两年的连续监测，揭示该河道弯道的岸滩侵蚀后退主要发生在冬季［10］。该研究表明冬季河道水沙运动过程对岸滩演变具有重要影响，但没有从机理上解释引起岸滩冬季侵蚀的原因。Zabilansky等在1998—1999年冬季对美国中部密苏里河的7个代表性断面开展现场观测，发现这些断面均存在显著的流凌刮擦侵蚀［11］。该研究揭示了河冰运动与冬季岸滩侵蚀后退的联系，解释了冬季河道岸滩侵蚀频发的诱因。他们的研究主要侧重于现场观察分析，并没有从机理上解释河冰运动对岸滩侵蚀的促进机制。Beltaos 等在加拿大阿萨巴斯卡河下游的冰塞观察中也发现冰塞体的释放能引起巨大的泥沙输运，并伴随着严重的河岸刮擦侵蚀［12］。这些研究均强调河冰运动对岸滩侵蚀的影响，但缺少对河冰与岸滩作用力的分析，也没有分析岸滩崩塌破坏机制。

关于河岸崩塌破坏的研究大部分集中在无河冰影响的工况。Simon等将土力学坡体稳定性分析中的静力平衡理论应用到河岸稳定分析，建立堤脚冲刷和岸坡崩塌调整模型［13-14］。该模型能成功模拟河道水流冲刷引起的堤脚淘刷和岸坡失稳变化过程，为堤岸稳定性分析提供了有效的数值模拟技术。Duan 和Julien 采用岸坡平行后退假设，建立耦合岸滩演变和水沙运动的平面二维数学模型，并成功用于长周期弯道演变机理研究［15］。该研究假设初始河岸坡度处于临界稳定坡度，但这种假设并不符合实际河道复杂的岸坡条件，不能准确模拟河岸演变过程［16］。Abderrezzak 等采用统一泥沙休止角判别岸坡稳定与否，在满足泥沙质量守恒的基础上调整岸坡坡面，并成功模拟实验条件下的河岸演变过程［17］，但统一的泥沙休止角并不能反映不同土体性质和含水量的岸坡稳定性［18］。Zech等采用两个不同的泥沙休止角分别代表水上和水下岸坡的临界失稳坡度，并采用较小的休止角重新分布崩岸土体［19］。该方法适合自然环境下非均质岸坡的稳定分析，成功模拟了实验条件下溃坝引起的岸坡崩塌侵蚀过程。这些研究为岸滩侵蚀和崩塌破坏提供了有效的技术手段，但没有涉及河冰运动对岸坡稳定性的影响分析。

综上所述，目前的研究缺少对河冰运动与岸滩侵蚀相互作用的分析。本文提出一种北方河流二维水冰沙耦合数学模型，旨在研究冬季北方河流复杂的水冰沙相互作用机制。研究采用耦合水动力过程、河冰运动过程、泥沙运动、河床演变和岸滩侵蚀的方法，模拟了实验条件下溃坝引起的岸坡崩塌侵蚀过程、冰盖下的岸滩侵蚀及冰塞冰坝形成与释放引起的岸滩崩塌和河床变形，有利于揭示河冰运动与岸滩侵蚀间的耦合作用机制。

2 水冰沙耦合数学模型

以Hung Tao Shen河冰研究团队开发的二维河冰水沙动力学模型为基础［20-21］，潘佳佳等建立了一种北方河流二维水冰沙耦合数学模型［22-23］。该模型耦合了二维水沙数值模块、河冰动力学模块和岸滩侵蚀模块，能模拟冰塞冰坝等极端工况下的水位流量波动、河冰聚集和释放、流凌对岸滩和河床的刮擦侵蚀、泥沙输移和岸滩的崩塌侵蚀过程［24］。模型的框架示意见图1。在给定地形、水位、流量和气温等初始条件和边界条件下，水冰沙耦合模型先采用具有迎风特性的Petrov-Galerkin（SUPG）型有限元法计算三角形非结构网格上的水位、流速、流量、水温、泥沙运动和河床变形，然后将计算结果传递给河冰动力学模块［25-30］；利用无网格的光滑粒子法（SPH）计算河冰分布、冰速、冰厚、面密度和冰浓度等信息，将计算结果传递给岸滩侵蚀模块［23］；采用双泥沙休止角法判断不同含水层岸坡的稳定与否，在质量守恒的基础上调整河岸崩塌后的坡面及失稳土体堤脚的再分布；最后将调整后的岸滩边界反馈给河冰动力学模块和二维水沙数值模块，校正地形变化后的水位流量、河冰运动、泥沙输运和河道变形。基于以上3个模块间的信息传递和数据反馈实现水冰沙的耦合模拟。

图1 二维水冰沙耦合数学模型框架

二维水冰沙耦合数学模型旨在解决河流全季节的水沙运动、河床及河岸演变问题。其中，二维水沙数值模块负责计算河冰影响下的水位流量波动、泥沙输运及河床冲淤变化，河冰动力学模块主要模拟水流驱动下的河冰运动、堆积、冲蚀及释放过程，而岸滩侵蚀模块提供流凌刮擦下的河道边界变化。该模型适用于以下研究：（1）复杂河道边界和地形突变下的急缓流交替过程；（2）降雨降雪、太阳辐射和不同风场气温影响下的水温变化过程；（3）冬季低气温影响下的水内冰、锚冰、岸冰和浮冰的生长、消融及变化规律；（4）冰盖下流冰花的侵蚀、输移和沉降过程；（5）开河和封河过程中的冰塞冰坝产生、发展和释放机制；（6）河流冰盖的热力增长消融和武开河过程；（7）河冰影响下的泥沙运动和河床侵蚀规律；（8）流凌刮擦和冰塞冰坝影响下的岸滩崩塌侵蚀及再平衡过程。具体的理论背景和基本方程参考北方河流水冰沙耦合研究［24］。

水冰沙耦合数学模型考虑了河冰的形成、发展、输移、堆积、释放及热力生长消融过程，能计算冬季河流水和冰的相互作用，也能研究夏季或无冰河流水、沙、河床与河岸间的相互作用机理，还能用于北方河流全季节的水冰沙耦合过程模拟，将水沙理论和河冰理论相融合［22］。

3 模型验证

二维水冰沙耦合数学模型［20-21］的河冰动力学模块已被多个实验算例验证，并成功用于黄河、加拿大皮斯河和美国尼亚拉加河等多个河流的河冰过程模拟［20，25，28］。本文进一步验证了二维水冰沙耦合数学模型的岸滩侵蚀模块，重点采用Zech等试验条件下溃坝引起的河岸侵蚀后退资料检验该模块的准确性与可靠性［23］。试验设置在一均匀顺直水槽，断面形态为半梯形断面，具体地形与试验布置见图2。河床为平坡，水槽长12 m，宽0.495 m。梯形断面底部宽0.15 m，初始边坡为50°，坡高为0.16 m。初始时刻闸门设置在距离上游边界5 m处，位置见图3。闸门上游河道和河岸初始时刻是不可冲刷的，闸门开启后可冲刷，水位高0.15 m；其下游是干河床，河道与河岸由1.8 mm的均匀沙组成，相对容重为2.615。该均匀沙在干燥状态下泥沙休止角为82°，大于初始边坡角度；其水下泥沙修止角为30°，小于初始边坡角度。因此岸坡在无水条件下是稳定的，在水位淹没以下部分是不稳定的。该试验展示了溃坝水流淹没引起的岸滩崩塌侵蚀过程，被用于河岸侵蚀模型的验证，具有较好的代表性。

图2 溃坝试验布置（改自Zech等［19］）

在试验中，闸门设置在x=0m处，闸门下游0.5 m 和1.5 m 分别设有断面地形高程记录仪，用于提供不同时刻断面侵蚀后退的试验数据。二维河冰动力学模型采用非结构的三角形网格，以适应自然河流复杂的地形边界变化。沿流线方向网格精度为0.05 m，横断面方向网格精度为0.016～0.1 m。网格精度在关心的坡面区域为0.016 m，在主河道区域相对较粗。计算的时间步长为0.001 s。

图3展现不同时刻模拟的溃坝水位过程。当闸门在t=0s突然开启后，水跌引起的溃坝水波快速向下游传播，在t=1s时经过第一个测站0.5 m处；在t=3s时，水流头部到达第二个测站1.5 m处。5 s后水流传播到闸后3 m，但并没有遇到下游边界。因此下游边界的影响可以忽略不计。

图3 不同时刻模拟的溃坝水位过程

图4对比了本文模型、Zech 等［19］和Swartenbroekx 等［31］计算的断面侵蚀崩塌过程与实测资料。在t=3s时，Zech等模拟的坡面侵蚀过大，而堤脚附近淤积偏少；Swartenbroekx等模拟的坡面侵蚀也偏大，而堤脚淤积偏小。本文所采用模型计算的0.5 m和1.5 m两个横断面地形高程与实测资料吻合良好，模型计算精度优于其他两个模型。在t=5s时，Zech 等模型计算的两断面坡面崩塌面积大于实测值；Swartenbroekx 等模拟的坡面侵蚀面积较试验值偏小，堤脚淤积量也偏小。水冰沙耦合数学模型模拟的断面形态与试验资料最为吻合。对比分析可知，本文提出的模型能准确模拟溃坝水流引起的岸坡侵蚀和崩塌过程，而Zech等和Swartenbroekx等模拟的结果存在较大误差。这主要是因为Zech等采用一维水沙数值模型，而断面形态仅通过5个点来捕捉坡面侵蚀和堤脚淤积，横断面的网格精度较低。Swartenbroekx 等采用平面二维数学模型，但其输沙模型采用平衡输沙模式，不能准确模拟溃坝水流引起的急剧坡面侵蚀和崩塌破坏过程。本文采用非平衡的水沙运动模块、河冰动力学模块及岸滩崩塌侵蚀模块相耦合的数学模型，能更准确地揭示岸滩侵蚀、失稳坡面变化和坍塌土体再分布规律。

图4 不同时刻和位置下模拟的岸坡侵蚀过程与实验资料的比较

4 算例应用

4.1 恒定均匀流模拟在溃坝水槽实验尺度下验证二维水冰沙耦合数学模型后，该模型被进一步应用于实际尺度冰盖冰塞影响下的河床冲淤和岸滩侵蚀研究。概化的数值计算区域为顺直、均匀梯形河道，代表性横断面形态及正常水位见图5。河底宽12 m，计算区域宽34 m，河岸坡度为1∶4.33，坡顶高度为2.2 m。河道长为400 m，河床底坡为0.0005。河岸与河道由粒径为0.5 mm 的均匀沙组成。计算区域的非结构三角形网格，共计节点5025个，三角形单元9600个。沿流线方向网格精度为2 m，横断面方向网格精度在0.9～2 m 间，其中坡面处网格精度较高。网格沿河道中心线y=17m对称，以保证非结构网格计算的准确性。数学模型计算的时间步长为0.05 s，模型所采用的主要参数见表1。上游边界设为恒定的流量30m3/s，下游边界条件设为正常水深1.6 m。在忽略泥沙运动和河冰运动时，模拟的恒定均匀流水深和流速分布如图6所示。模拟的水深和流速沿轴线对称，且分布均匀，验证了该模型在均匀对称河道水流计算中的可靠性。断面平均流速为0.98m/s，弗劳德数为0.29。该恒定均匀流模拟为以下河冰运动和岸滩侵蚀模拟提供了可靠的初始条件。

表1 河流水冰沙耦合模型关键参数

图5 均匀梯形断面形态及水位示意

图6 恒定均匀流下模拟的水深和流速分布

4.2 冰盖对河道冲淤的影响在图6计算的恒定均匀流下，分别设置有无冰盖的计算工况，上游泥沙边界采用循环边界，即下游出口流出多少泥沙，上游边界就提供多少泥沙。循环泥沙边界适用于长河道的水沙模拟。图7为明渠和冰盖覆盖两种工况下100 h 和200 h 模拟的水位和断面地形。模拟200 h后梯形河道达到冲淤平衡状态，冰盖条件下稳定的水位高于明流工况。这主要是因为冰盖增加了河道边界的湿周和总体阻力，河道断面平均流速下降，因此考虑冰盖时的河道水位更高。4个不同断面均显示冰盖下水位更高，水面下岸坡的冲刷起点更高，但冰盖下岸坡总体冲刷量小于明流工况，相应主河道的泥沙淤积也偏小。这主要是因为冰盖下的流速比明流工况低，岸坡和床面的水流拖曳力低于明流工况，因此岸坡的冲刷和主河道的淤积小于明流条件。对比100 h和200 h的模拟结果，下游河道比上游更快达到冲淤平衡。计算结果显示二维水冰沙耦合数学模型能模拟河道中冰盖对岸坡冲刷和河道淤积的影响。

图7 明流和冰盖影响下不同断面模拟的水位及河床冲淤对比

4.3 冰塞形成和释放过程为了进一步验证二维水冰沙耦合数学模型在冰塞冰坝等极端工况下的可靠性，研究采用第3节工况下的梯形断面河道及河岸泥沙条件，模拟了冲积河流冰塞冰坝的形成和释放过程。在河道中间设置拦冰栅并假设冰塞形成过程中河床为可冲刷而河岸固定，在冰塞释放过程中河床与河岸均可冲刷。图8显示了稳定冰塞形成后浮冰颗粒的分布、河道地形的冲淤变化、纵剖面的水位、上下冰面及河床高程分布。图中，黑色方块示意浮冰，地形冲淤变化中蓝色示意冲刷，红色示意淤积，两条白线示意左右岸堤脚位置。模拟结果显示卡冰处冰体堆积严重，最大冰厚超过0.12 m，冰塞引起上游水位壅高，上下水位差达0.04 m。此外，冰塞体前端冰厚最大，并导致断面过流面积减小，局部流速增大，河道主槽出现显著冲刷。模拟结果与Beltaos等在加拿大阿萨巴斯卡河观测的冲积河流冰塞形成过程相符合［32］，证实了模型对冰塞形成中水冰沙耦合模拟的可行性。

以图8模拟的冰塞体为初始条件，在t=0 时刻撤去拦冰栅，进一步采用二维水冰沙耦合数学模型模拟冰塞释放过程中的河冰运动、河道冲淤变化、岸坡崩塌和再平衡、水位流量波动过程。模拟的不同时刻水位、地形和河冰分布结果见图9。计算2 s后，储存在卡冰点上游的浮冰和水体以类似溃坝水体的洪水波向下游传播。下游岸坡在淹没和水流冲刷下发生河岸侵蚀、崩塌和堤脚淤积，而上游河道流凌及水流也导致岸坡上部冲刷坍塌及堤脚淤积。随着冰塞的释放，冰塞体前端水位及冰厚减小，并向下游平滑过渡。在3至6 s，上游储存的浮冰和水体继续向下游释放，流凌刮擦引起的岸滩崩塌更为显著，坍塌的土体堆积于堤脚附近，部分泥沙被输移到下游。计算的12 s，冰塞释放的浮冰流出计算区域，上游壅水部分回落，水面附近的岸坡出现显著侵蚀和崩塌，而两侧堤脚出现大量淤积。其中上游的岸坡冲刷和堤脚淤积大于下游河道。这是合理的，因为上游河道壅水高于下游，冰塞体刮擦影响区域也较下游河道大。

图8 冰塞形成引起的水位壅高及河床冲刷模拟结果

图9 不同时刻模拟的冰塞释放过程及水位、地形和河冰分布

图10显示了上游至下游4个典型断面0、3、6、9和12 s计算的断面地形冲淤变化过程。4个断面均显示水面下存在岸坡冲刷侵蚀，水面以上的岸坡也出现崩塌侵蚀，大量失稳土体落淤在堤脚附近。从上游到下游，冲刷和崩岸的土体面积逐渐减小，而堤脚淤积的土体面积相应也减少。这主要是因为水位从上游到下游逐渐减小，河冰大量堆积在上游河道，流凌刮擦和水流侵蚀范围从上游向下游递减。此外，水面以上岸坡再平衡的坡角大于水面以下岸坡的平衡坡角。计算结果显示水冰沙耦合模型能准确模拟水面上下不同含水层的稳定坡面。图8—10 证实该模型能有效模拟极端冰塞下的岸滩刮擦侵蚀、岸坡崩塌、失稳土体堤脚落淤和冰塞释放过程，能进一步应用于其它水冰沙问题研究。

图10 冰塞释放引起的岸滩侵蚀和断面地形变化

5 结论

本文初步建立二维水冰沙耦合数学模型，由水沙数值模块、河冰动力学模块和岸滩侵蚀模块3部分组成。模型采用非结构的有限元法计算水沙运动和河床冲淤变化，能适应自然河流复杂的河道边界；采用无网格的光滑粒子法计算河冰运动、堆积和释放过程，能模拟极端冰塞冰坝的形成、发展和释放过程；采用双泥沙休止角法计算岸滩的侵蚀、崩塌和再平衡过程，能揭示岸坡的周期性侵蚀规律。研究采用溃坝水流引起的岸滩侵蚀和崩塌实验验证了模型的准确性和可靠性，并将模型应用于冰盖和冰塞形成与释放引起的河床冲淤与岸滩崩塌侵蚀分析。模拟结果显示二维水冰沙耦合数学模型能揭示流凌对岸滩的刮擦侵蚀、河岸崩塌及失稳土体再分布规律，为极端冰塞冰坝下的水冰沙相互作用研究提供了新工具。二维水冰沙耦合模型创新性地融合河冰理论和水沙理论，同时满足无冰和有冰河流水沙和河冰模拟需求，为河流全季节的水冰沙耦合问题研究提供了有力支撑。

北方河流冬季水流运动、河冰输运、泥沙输移、河床变化与岸滩侵蚀存在复杂的水冰沙相互耦合作用，本文在数值模拟和理论计算方面开展了详细研究，还需通过原型观测和物理模型实验等多种手段深入研究水冰沙耦合作用机理。关于北方河流河冰与河岸的相互作用及泥沙运动的实测资料和理论分析欠缺，亟需开展更多系统水冰沙耦合作用研究。