中国北方季节性冰冻河流岸滩崩塌数值模拟——以松花江为例

假冬冬,杨 俊,陈长英,张幸农,应 强

(南京水利科学研究院水文水资源与水利工程科学国家重点实验室,江苏 南京 210029)

岸滩侧蚀崩塌现象广泛分布于世界各大江河之中[1-3],是一种危害较大的自然现象,与水沙动力等因素共同引起复杂河床冲淤调整[4]。对于中国东北地区的季节性冰冻河流(季冻区河流)而言,河道岸滩、岛屿汊道不稳定,冻融作用强、崩岸频发[5-6],直接引起农田土地流失。此外,季冻区河流受水动力-冻胀/冻融耦合作用,其岸滩崩塌与一般河流存在明显差异,机理更为复杂。因此,开展季冻区河流岸滩崩塌机理及其数值模拟研究,可为江河治理与保护提供科技支撑,具有重要的理论和实践意义。

岸滩崩塌是指在近岸水沙与河床边界的相互作用下,河岸受到各种因素影响而发生部分土体崩塌的现象[7]。随着对崩岸现象的深入认识,国内外学者在崩岸机理与力学模式方面开展了具有重要价值的研究工作,建立了不同类型土质河岸的岸滩崩塌描述模式[8-11],为河道崩岸的数值模拟提供了力学基础。基于此,不少研究者将岸滩崩塌力学模式与传统水沙输移模型相耦合,开展了岸滩侵蚀过程的二维数值模拟和三维数值模拟研究,有效复演了岸滩崩塌及河道演变过程[12-15],并在长江、黄河等河流岸滩失稳过程中得到成功应用。此外,传统水沙数学模型在松花江等季冻区河流河道治理中也得到广泛应用,例如:宗原等[16]对黑龙江雪水温河段护岸的防护效果进行了研究,护岸的修建可有效控制河势稳定，防止因河岸崩塌引起的国土流失；陆永军[17]采用二维泥沙数学模型,模拟了松花江干流三姓浅滩航道整治一期工程丁坝群作用下的回流边线、汊道分流比以及推移质运动引起的河床变形,但此类模型并未考虑季冻区河流岸滩崩塌机理。从当前岸滩崩塌机理及其数值模拟的研究进展来看,主要以水动力冲刷下岸滩崩塌机理及其与河床冲淤交互作用为主,尚未深入考虑水动力-冻胀/冻融耦合作用,但对于季冻区河流而言,受气候影响则存在显著冻胀/冻融特点,且对岸滩力学强度和稳定性具有重要影响[18-20]。除此之外,对于松花江等季冻区河流而言,冻融作用强,岸滩崩塌现象常见,但受模拟技术的限制,其河道演变模拟尚未考虑岸滩崩塌过程,亦未涉及冻胀/冻融作用对岸滩崩塌的影响。

本文针对季冻区河流岸滩崩塌特点,以松花江干流大顶子山航电枢纽下游近坝段为例,在三维水沙数学模型研究成果[14-15]的基础上,通过建立水动力-冻胀/冻融耦合作用岸滩崩塌力学描述模式,提出季冻区河流岸滩崩塌及河床变形的数值模拟方法,采用实测数据对数学模型进行验证,在此基础上分析岸滩崩塌对河床冲淤的影响。研究成果可为季冻区河流岸滩崩塌及河道演变的深入研究提供技术支撑。

1 水动力-冻胀/冻融耦合作用岸滩崩塌机理

对于季冻区河流而言,岸滩侵蚀及其演变过程除了受水沙条件影响外,因气候影响还存在显著冻胀/冻融特点,且对岸滩稳定性造成重要影响,机理更为复杂,传统河流动力学理论已不能完全适用。每年封冻初期,河岸土体内的水体在负温度条件下结冰,土在冻结时易发生体积膨胀,即发生冻胀现象。冻胀作用破坏了河岸土层的结构,土壤结构变松,经常使得岸滩土层出现裂缝,甚至坍塌。而到次年春天,当上升到正温度时,因冻结土层的冰融化,作为土体临时结构黏结力的冰晶体消失,土体结构变松,力学强度降低,对于处于极限平衡状态的岸坡而言,冻融作用会进一步使河岸边坡处于不稳定状态,最终导致河岸坍塌。

基于上述分析可知,冻胀作用通过冻胀力会使岸坡开裂,而冻融过程则会使岸滩土体力学强度降低,两者均会影响岸坡的稳定。因此,对于季冻区河流岸滩崩塌而言,既要考虑水动力对岸滩的冲刷侵蚀,又要考虑冻胀/冻融对岸滩稳定的影响机理。因此,需建立坡脚水力冲刷计算模式以及冻胀/冻融耦合作用下岸滩失稳模式。

1.1 坡脚水力冲刷计算模式

水流淘刷作用下,坡脚侧蚀冲刷过程可由笔者提出的岸滩侧蚀速率公式计算[11]:

(1)

式中:ωb为岸滩侧蚀速率；λ为侧蚀系数,由实测资料进行率定；γ为水体重度；γb为岸坡土体重度；u为近岸处水流流速；uc为泥沙起动流速。

1.2 岸滩失稳模式

1.2.1 上部土体裂缝深度计算

在水力冲刷与冻胀作用下,坡顶开裂,会形成一定深度的裂缝,并影响岸坡稳定。根据土力学朗肯理论,考虑冻胀作用时,岸滩坡顶裂缝深度为

(2)

式中:Ht为考虑冻胀作用的坡顶裂缝深度；pd为冻胀力；c为土体黏聚力；Kz为朗肯主动土压力系数。

冻胀力与土体性质(黏聚力)密切相关,设ξ=pd/(2c),则式(2)可简化为

(3)

分析认为,冻胀作用会使岸滩坡顶开裂深度进一步增加,从而降低岸滩的稳定性。

1.2.2 岸滩崩塌力学模式

水动力-冻胀/冻融耦合作用下岸滩崩塌过程,同样基于岸滩受力分析,可推导获得岸滩崩塌力学模式。从三维受力分析角度,考虑沿纵向相邻土体对其稳定性的影响,笔者借鉴Osman和Thorne[21]的思路,推导建立了岸滩崩塌力学模式[15]。与此类似,考虑冻胀对土体上部开裂深度影响以及冻融对土体力学强度影响,得到水动力-冻胀/冻融耦合作用下岸滩崩塌力学模式。

(1)初次崩塌。在已知季冻区河流初始岸滩高度(H0)、初始岸坡(i)的情况下,根据水沙数学模型可计算得到坡脚横向冲刷值为ΔB,河床冲深为ΔZ，并根据垂向、横向变形计算结果,修正岸滩形态。图1为发生冲刷后初次崩塌的示意图。初始岸滩地形由ABCD连线表示，经水流冲刷及崩岸后的岸滩地形变为A′B′C′D。FS为土体滑动力,FR为土体抗滑力,W为岸滩单位宽度坍塌土体的重量,WB为坍塌土体的宽度,β为发生初次崩塌时,滑动面与水平面的夹角,可由式(4)计算:

图1 岸滩初次崩塌Fig.1 Sketch of the initial bank failure

(4)

式中:H为冲刷后的岸滩高度；H′为岸坡线转折点以上的岸滩高度；k=Ht/H；i为初始岸坡；φ为岸滩土体的内摩擦角。

对于纵向长度为Δl的岸滩土体,假设FS=FR时,岸滩处于崩塌的临界状态,由此可求得岸滩崩塌的临界条件为:

(5)

式中:Δl为岸滩崩塌的纵向长度,pi-1、pi+1分别为上、下游侧面所受黏聚力系数。

在计算得到临界值(H/H′)cr和实际值(H/H′)m的情况下,按如下方法可判断岸滩是否会崩塌:当(H/H′)m<(H/H′)cr时,岸滩稳定,不是发生岸滩崩塌的临界条件,进入下一个时段的水沙计算；当(H/H′)m≥(H/H′)cr时,岸滩失稳,利用几何形态关系,可计算坍塌土体尺寸。

(2)二次崩塌。岸滩发生前述初次崩塌后,并假定之后的岸滩崩塌方式为平行后退,即岸滩崩塌的破坏角仍为β,崩塌示意图见图2。

图2 岸滩二次崩塌Fig.2 Sketch of the subsequent bank failure

此时,假设FS=FR时,可获得岸滩二次崩塌临界条件为

(6)

与初次岸滩崩塌判别类似,在计算得到临界值(H/H′)cr和实际值(H/H′)m的情况下,即可判断岸滩是否发生崩塌。

(3)冻融作用引起的岸滩土体力学性能变化。土体冻融过程中,因冻结土层的冰融化,作为土体临时结构黏结力的冰晶体消失,土体结构变松,强度降低。其中，抗剪强度是对岸滩稳定性具有重要影响的力学性质指标,通常采用黏聚力和内摩擦角来表示。根据已有研究表明,冻融循环作用下土体黏聚力显著减小,但内摩擦角的变化相对不明显[18-19]。根据文献[19]的研究,多根样品曲线试验结果表明:黏聚力随着冻融循环次数的增加而逐渐减小,单次冻融后黏聚力的减幅范围约为15%～40%；而内摩擦角无规律可循呈上下波动的状态，单次冻融后内摩擦角有增有减,变幅不大。根据文献[20]的研究,单根样品曲线试验结果表明:单次冻融后黏聚力的减幅约为19%；单次冻融后内摩擦角的增幅约为14%。因此,冻融作用对岸滩土体抗剪强度和稳定具有重要作用,本次模拟研究中将考虑冻融影响,对抗剪强度进行相应调整。

松花江干流属于典型的季节性冰冻河流,每年4月至11月一般为畅流期,造床作用相对较大,12月至次年3月一般为封冻期,河床冲淤相对较小。针对这一特点,本研究主要对畅流期进行模拟,计算过程中对其进行概化处理:4月考虑冻融影响,即抗剪强度相应调整；11月考虑冻胀作用,即式(3)中考虑冻胀作用,坡顶开裂深度增加,从而降低岸滩的稳定性；其余时段模拟则考虑常规岸滩稳定计算。

2 岸滩崩塌与河床冲淤数值模拟技术

岸滩崩塌与河床冲淤的模拟,一方面需要对传统水沙输移和河床垂向冲淤过程进行计算,另一方面还需对岸滩崩塌过程进行模拟,并实现两者的有机联系。本节将对季冻区河流岸滩崩塌与河床冲淤的模拟技术进行介绍。

2.1 局部网格可动技术

与常规河流模拟类似,在应用数学模型对季冻区河流岸滩崩塌及河床冲淤过程进行模拟时,崩岸的发生会改变数值模型的地形和边界,而岸滩崩塌宽度不一定与崩岸处的模型网格宽度一致,计算网格很难准确跟踪崩岸后的地形和边界。本文采用一种基于非正交网格的局部网格可动技术对岸滩崩塌过程进行岸坡边界的跟踪,在整个大计算域内生成计算网格。在模拟过程中,仅对岸滩崩塌附近的网格进行移动,实时准确地跟踪新的岸坡边界,其余计算网格位置保持不变。这样做既可较准确地实时跟踪岸滩崩塌后的岸坡边界,又无需重新生成整个计算域内的模型网格,可弥补传统定网格和动网格在这方面存在的不足之处。具体过程可详见文献[14]。

2.2 数值模拟流程

本文模拟过程中采用的三维水流泥沙模型为非均匀非平衡泥沙输移模型,模型基于非正交曲线网格系统,采用控制体积法进行方程离散与求解[14-15],该模型已在实际天然河流中得到了较好的验证和应用。岸滩崩塌模拟过程中,首先给定各变量初值,采用三维水沙模型计算河道内流场、水位等水动力信息；依据水动力模拟结果,计算悬移质含沙量及推移质输沙率；随后计算河床冲淤变形,并更新地形；在进行冻融/冻胀作用判断与计算的基础上,根据岸滩崩塌机理与力学模式,计算和判断岸滩稳定状态,若岸滩失稳,则修改岸滩及近岸地形以及河床级配信息；依次重复上述计算过程,直至完成设置的计算时间。

3 典型季冻区河流岸滩崩塌数值模拟

以松花江干流大顶子山航电枢纽下游近坝段为例,进行典型季冻区河流岸滩崩塌数值模拟分析。松花江干流属于典型的平原冲积河流,全长约940 km,区间集水面积18.64万km2,两岸除部分区段为丘陵岗地外,基本为冲积平原,向东流至同江附近汇入黑龙江。大顶子山航电枢纽工程位于松花江干流中游、哈尔滨市区下游约63 km处,是松花江航道梯级开发总体规划中8座枢纽之一,工程于2004年9月开工建设、2008年末交工使用。枢纽的运行以及上游来沙减小,使得枢纽下游近坝段出现较为明显的冲刷,深槽摆动,洲滩崩退现象时有发生。本节将依据前文模拟技术,建立大顶子山航电枢纽下游近坝段河道演变动力学模型,并采用实测资料对模型进行验证,分析岸滩崩塌对河道演变计算结果的影响。

3.1 计算范围及计算条件

计算河段上起大顶子山航电枢纽,下至富江岛汇流口下游顺直段,模拟河段全长约22 km。基于2008年10月和2014年10月实测河床地形构建数学模型,模型采用贴体曲线非正交网格与复杂岸线贴合良好,模型平面计算网格为880×280,平均网格尺度沿水流方向约为25 m、沿断面方向约为10 m,并在局部区域对计算网格进行适当加密,使网格能够较好地反映复杂河道边界,模型垂向网格分为13层。模型平面计算网格划分及河势分别见图3、图4。模拟时段内模型进口流量过程如图5所示,出口水位由率定的水位—流量关系确定(图6)。松花江干流属于少沙河流,主要由推移质造床为主,且计算河段位于大顶子山航电枢纽下游近坝段,枢纽运行初期模型进口不考虑上游悬沙。

图3 计算网格划分Fig.3 Computational meshes

图4 研究河段河势Fig.4 River regime of the study area

图5 模型进口流量过程Fig.5 Flow discharge processes of model inlet

图6 计算水位—流量关系与实测水位—流量关系对比Fig.6 Calibration of relationship between water level and discharge

3.2 模型计算分析与讨论

(1)水位—流量关系验证。为使模型能够准确复演天然河道水流运动状态,首先采用模拟河段进口(大顶子山航电枢纽坝址)及下游(模型出口处)的水位—流量关系进行计算,并与天然水位—流量关系进行比较。各级流量下模型计算水位—流量关系与天然水位—流量关系对比情况见图6,由图可见,模型计算水位—流量关系与天然情况基本吻合,各级水情下水位差值最大一般不超过0.10 m,说明模型能够较好地复演河道阻力情况,率定的河段综合糙率一般为0.022～0.035,主槽阻力相对较小,滩地阻力相对较大(综合糙率一般在0.030～0.050左右)。

(2)典型流场分布特征。本河段典型流量下(多年平均流量Q=1 200 m3/s)流场分布特征如图7所示。由图可见,中水期河道内水流流速一般在1 m/s以内,顺直段流速较为平顺,弯曲段三维流场特征较为明显,表层、中层流速大于底层流速,表层流速指向弯道凹岸、底层流速指向弯道凸岸,局部最大夹角可达40°左右。江中的富江岛将该处河道分为左右两汊,其中左汊为支汊、右汊为主汊,两汊分流之比约为1∶4；两汊河道均较为弯曲,凹岸侧流速较凸岸侧大,弯顶处尤为明显,易造成河道凹岸冲刷。从三维流场中截取典型断面(CS1、CS2)流速分布(由下游往上看)如图8所示,断面位置见图7(b)；从图中可看出,弯道断面环流特征较为明显,表层流速流向凹岸、底层流速流向凸岸。

图7 流速分布特征Fig.7 Computed velocity distribution

图8 典型断面流速分布特征Fig.8 Computed cross-sectional velocity distribution

(3)典型断面冲淤验证分析。根据实测河道地形资料,对研究河段的河床冲淤变化进行验证。收集了计算河段2008年10月实测地形资料,作为模型计算初始地形,根据2014年10月实测的典型大断面资料对模型进行验证。根据现场采样分析,岸滩土体中值粒径约为0.09 mm,级配曲线如图9所示,黏聚力以及重度分别取值为14.0 kN/m2和18.0 kN/m3,内摩擦角取值23°,模拟过程中采用此值作为基准值。冻融期对抗剪强度基准值进行相应调整,本文首先根据文献[19]的研究进行调整,即单次冻融后黏聚力减幅为30%(取均值),内摩擦角保持不变。

图9 岸滩土体级配Fig.9 Sediment gradation of riverbanks

大顶子山航电枢纽运行初期,依据收集到的资料,下游近坝段典型断面冲淤变化验证见图10,断面位置见图4。由图可见,总体而言2008年10月至2014年10月各断面以冲刷为主,近坝处断面冲刷尤为明显,断面平均冲深在2 m左右,最大冲刷在4 m以上,下游弯道处断面弯道凹岸冲刷明显,岸线崩退,最大崩退距离可达150 m左右,弯道凸岸则出现淤积。具体而言,从各典型断面的冲淤变化来看,DM1位于枢纽下游约1.5 km处,枢纽运行初期(2008—2014年)河床在整体冲刷下切的同时,右侧深槽相邻河岸崩退距离约100 m；DM2位于枢纽下游约3.0 km 处,与DM1类似,河床整体冲刷下切约2.0 m,右侧河岸后退距离约150 m；DM3位于枢纽下游约6.0 km处,河床表现为左侧冲刷,右侧淤积,总体上仍呈现为冲刷,平均冲刷1.3 m左右；DM4位于枢纽下游约10.0 km处,位于弯道进口段,受水流取直影响,右侧岸坡不断崩退,崩退距离约200 m,河道展宽明显；DM5位于枢纽下游约12.0 km处,位于弯道段,受弯道环流和水流顶冲影响,凹岸侧河床冲刷、岸线后退约140 m,凸岸侧河床则出现淤积；DM6位于富江岛右(主)汊弯道处,凸岸出现淤积,凹岸则冲刷明显；DM7位于富江岛两汊汇流处下游,河道中部冲刷,两侧有所淤积。由上述各断面冲淤验证结果可看出,考虑岸滩崩塌的模拟结果与实测结果冲淤定性基本一致,冲淤幅度基本相当,数学模型可较好地反映河道的冲淤变化、岸滩崩退以及河道展宽过程。

图10 典型断面冲淤对比分析Fig.10 Simulated and measured results of deposition and erosion distribution at typical cross-sections

(4)岸滩崩塌对断面冲淤影响分析与讨论。为分析岸滩崩塌对断面冲淤的影响,对考虑与不考虑岸滩崩塌的模拟结果进行了对比,见图10。从结果对比来看,两者存在较明显的差异,考虑岸滩崩塌的模拟结果与实际情况相对吻合,由此也说明岸滩崩塌对于季冻区天然河道演变准确模拟的重要性。相对而言,顺直段两者差别较小,河床主要以冲刷下切为主,弯曲段受水流取直和顶冲作用,岸滩冲刷后退现象较明显,考虑岸滩崩塌后的模拟结果与实际观测较为一致,而未考虑岸滩崩塌的模拟结果则与实际情况存在较大差异。

为进一步分析冻融期抗剪强度调整对岸滩崩退的影响,除了前述根据文献[19]的研究进行调整外,即单次冻融后黏聚力减幅为30%(取均值),内摩擦角保持不变；还根据文献[20]的研究,对岸滩崩退进行模拟,即单次冻融后黏聚力减幅为19%,内摩擦角增幅为14%。前者简称工况1,后者简称工况2。对于岸滩崩退不明显的断面而言,两工况的模拟结果基本一致；为方便比较两者模拟结果之间的差异,选取具有典型岸滩崩退的DM4、DM5进行对比分析,见图11。由图可见:两工况均有较明显的岸滩崩退过程,其差异并不十分明显,主要是工况1的岸滩崩退距离较工况2略大,最大差异在15 m左右。这可能是因为冻融期,河道流量和流速相对较小,河床冲刷幅度亦相对小,岸滩形态主要由前期河床变形所主导,只有当前期岸滩已处于或接近于临界崩塌状态时,此时抗剪强度降低才易引发岸滩崩塌。

图11 冻融抗剪强度变化对典型断面岸滩崩塌影响Fig.11 Effects of freeze-thaw on bank erosion at typical cross-sections

同时也应指出,尽管考虑岸滩崩塌后的模拟结果与实测结果相对吻合,但具体定量上的冲淤分布还存在一定差异,一方面是因为季冻区河流岸滩崩塌机理极为复杂,涉及水、土、气候等方面因素,不仅存在平行崩岸类型的崩退,也存在淘刷落崩的类型,而本文只考虑了平行崩退模式,且实际土质条件与模拟设置的土体抗冲性存在差异,这会对模拟结果产生影响,其精细描述模式及其模拟方法还有待进一步深入研究；另一方面,在春季融化时期,北方河流普遍存在冰凌漂浮现象,虽然持续时间不长,但冰凌漂浮对岸滩稳定性也存在较明显的影响,本文模型尚未对此进行考虑；此外,季冻区河流演变除了畅流期的主要冲淤过程外,封冻期亦存在一定幅度的冲淤变化,且涉及明渠流与有压流的相互转换,本文仅考虑了畅流期冲淤过程,并未对封冻期冲淤过程进行模拟,这同样也会对计算结果造成一定影响,这一复杂过程的精细表达与模拟亦有待进一步深入研究。

4 结 论

基于水动力-冻胀/冻融耦合作用岸滩崩塌机理,以松花江干流大顶子山航电枢纽下游近坝段为例,建立典型季冻区河流岸滩崩塌三维水沙动力学模型,开展岸滩崩塌与河床冲淤的数值模拟分析,主要研究结论为:

(1)季冻区河流岸滩崩塌除受水沙条件影响外,因气候影响还存在显著冻胀/冻融特点,且对岸滩稳定性造成明显影响。坡脚侧蚀主要由水力冲刷所控制,冻胀作用通过冻胀力会使岸滩开裂深度增加；冻融过程则会使岸坡土体力学强度降低,影响岸滩稳定；通过力学分析,初步建立水动力-冻胀/冻融耦合作用岸滩坍塌失稳描述模式,为季冻区河流岸滩崩塌模拟提供力学依据。

(2)将传统水沙输移和河床垂向冲淤计算与岸滩崩塌模拟相耦合,提出季冻区河流畅流期河道演变过程的数值模拟方法。岸滩崩塌对季冻区天然河流演变具有重要影响,考虑岸滩崩塌的模拟结果与实测结果冲淤定性基本一致,冲淤幅度基本相当,数学模型可较好地反映河床冲淤、岸滩崩退以及河道展宽现象。建立的季冻区河流岸滩崩塌三维数值模拟方法,为深入研究季冻区河流岸滩崩塌及其对河床演变的影响提供了一种有效技术手段。