窝崩是长江中下游常见的一种河道崩岸形式, 具有发展速度快、破坏力强等特点, 往往数小时即可塌进上百米的距离, 甚至破坏江堤, 严重威胁人民的生产生活, 是一种危害极大的自然灾害。2008年11月18日于长江南京河段栖霞区龙潭三官村友庄圩发生的三江口窝崩, 形成长约340 m、崩退230 m的窝塘, 导致长约200 m的主江堤水毁^[1]。因此, 在掌握窝崩机理的基础上, 深入研究窝崩形成过程的数值模拟方法, 可为崩岸预测与预警以及防护措施的确定提供技术支撑, 具有重要意义。

窝崩是水流冲蚀岸滩引发河道平面变形的一种具体表现形式, 其发展速度快且过程复杂, 现场实测资料非常有限, 目前多采用概化试验的方法对其形成过程进行研究^[2-3], 研究表明窝崩口门内的水流基本是脱离主流而在崩窝内形成的回流^[2], 窝塘内出现高强度的回流, 不仅对边壁形成淘刷, 使得土体失稳, 回流也不断地将滑落土体分散并搬运至窝塘外^[3], 随后回流逐渐减弱而形成新的泥沙淤积^[2]。传统冲积河流河床冲淤变化的数值模拟, 一般是针对垂向的冲淤变化, 并且依据冲积泥沙(即新淤土)运动规律基础理论^[4], 随着崩岸机理的深入认识^[5-10], 不少学者将传统水沙模型与崩岸力学模式相结合, 开展了类似条崩的岸滩侵蚀过程的平面二维数值模拟研究, 即通过二维水沙数学模型计算与河岸侧蚀崩塌力学模式相耦合, 成功地模拟了岸滩侵蚀及其相应的河床冲淤变化过程^[11-16]。为了模拟岸滩侵蚀的三维水沙输移过程, 文献[17-18]以长江中游荆江典型河段为例, 开展了三维水沙数模与二元结构河岸崩塌的耦合计算。而窝崩过程的水流结构、岸滩冲刷、泥沙输移均十分复杂, 回流环流等三维特征显著, 且回流区泥沙冲刷输移过程与常规明渠流存在较明显的差异^[19], 其数值模拟研究还有待进一步深入探索。

本文以2017年11月8日发生于长江下游扬中河段的指南村窝崩为例, 在三维水沙数学模型研究成果^[17-18]的基础上, 基于窝崩机理及其力学描述模式, 采用考虑回流区紊动影响的挟沙力公式模拟窝塘内泥沙输移, 建立典型窝崩的三维水沙动力学模拟方法, 并采用实测资料对模型进行验证分析, 在此基础上探讨窝塘内三维水流结构变化特征以及泥沙输移模式对窝塘形态模拟结果的影响。研究成果可为窝崩机理的深入认识以及窝崩治理提供科技支撑。

1   窝崩数值模拟技术

窝崩过程是长江下游常见的一种河床演变形式, 窝崩的模拟不仅需要对窝崩附近复杂水流结构进行计算, 还需对窝崩的机理及其力学描述模式加以考虑。

1.1   窝崩机理及力学模式

窝崩过程中, 一般伴有强烈的回流淘刷, 岸坡坡脚受到水流强烈的冲刷, 坡面土体冲蚀后, 其上部土体失去支撑, 平衡条件被破坏, 并随之产生坍塌, 坍塌的土体又被水流迅速分解搬运, 如此循环反复, 致使岸坡不断崩坍后退。因此窝崩模拟时, 需考虑三方面的过程:一是岸滩坡脚冲刷后退, 二是上部土体失稳坍塌, 三是水中土体的分解输移。

1.1.1   坡脚侧蚀冲刷计算

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

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

1.1.2   上部土体坍塌力学模式

窝崩过程中, 当下部坡脚淘刷到一定程度后(侧蚀宽度达到极限状态), 上部土体将会出现坍塌失稳。考虑纵向因素影响时, 上部挂空土层的临界挂空宽度ΔW_c为^[17]

式中：T₀、c、γ_b分别为上层土体的抗拉强度、凝聚力以及容重, 可通过实测得到；B、H_u分别为上层土体的宽度及高度；p_i－1、p_i+1分别为与其紧邻的上、下游侧面所受凝聚力的方向系数, 当其紧邻的断面上层土体判断发生坍塌时, 其取值为负, 否则取值为正, 绝对值范围为0~1。

根据下部坡脚冲刷后退距离ΔW以及上部土体临界挂空宽度ΔW_c即可判断上层土体的稳定与否:当ΔW < ΔW_c时, 岸滩上部土体不会受拉而发生坍塌；当ΔW=ΔW_c时, 岸滩上部土体处于坍塌临界状态；当ΔW>ΔW_c时, 岸滩上部土层将会受拉而发生坍塌。

1.1.3   坍塌土体分解输移模式

上部土体坍塌后, 一部分坍塌土体将堆积在坡脚河床处, 即按一定比例p_b均匀分布在近岸处^[18], 其余部分(1-p_b)则以源项形式转化为悬沙。

对于河道中的悬沙输移, 一般可采用如下形式的挟沙力公式^[21] :

式中: S_*为水流挟沙力；K、m为经验系数；g为重力加速度；R为水力半径；ω为泥沙沉速。该公式目前已被广泛接受, 并成功应用于长江等河道的实际工程中。

强烈紊动是窝崩快速发展的重要条件^[9], 对于回流区的悬沙输移过程, 可采用能量挟沙来推导^[19], 研究结果表明, 其紊动相对较强的区域, 水流输沙能力亦相对较强, 并存在如下关系^[19] :

式中: S_*,c为回流区水流挟沙力；η为回流区紊动强度与无回流区紊动强度之比。文献[19]通过试验研究得到, 丁坝后回流正流区和回流负流区的纵向最大紊动强度分别为无回流时的3.0倍和2.0倍, 由此说明紊动强烈回流区的水流输沙能力显著增强, 从而使得窝崩土体能够快速输移。

1.2   局部网格可动技术

在采用数值模型对窝崩过程进行模拟时, 窝崩的发生会改变模型的地形, 而崩塌宽度不一定正好与窝崩处的计算网格宽度相等, 这就使网格对边岸地准确拟合变得困难。本文采用基于非正交网格的局部网格可动技术对窝崩过程进行边岸的跟踪, 即在整个计算域内生成网格, 在模拟过程中, 仅对窝崩附近的网格进行移动, 使其能够准确地跟踪边岸位置, 同时其余网格位置不变。这样做既能较准确地拟合窝崩后的岸坡位置, 又无需重新生成整个计算域内的计算网格, 弥补了传统定网格和动网格在这方面存在的不足。具体过程详见文献[18]。

1.3   窝崩数值模拟流程

本文采用的三维水沙模型为非均匀非平衡泥沙输移模型, 模型采用非正交曲线网格建立, 采用控制体积法进行离散求解^[17-18]。窝崩模拟过程中, 首先给定各变量初值, 采用三维水沙模型计算该时刻流场、水位等水流信息；然后依据水流流场结果, 计算各粒径组悬沙含沙量及推移质输沙率；在此基础上, 计算河床冲淤变形, 更新地形；根据窝崩机理及力学模式, 计算坡脚侧蚀冲刷距离, 同时依据局部网格跟踪技术来识别窝崩后岸坡位置, 在此基础上, 依次判断上部土体稳定状态, 若失稳, 则修改岸坡及近岸地形；重复上述计算过程, 直至计算结束。

2   典型窝崩数值模拟

2.1   数学模型建立及验证

扬中河段位于长江下游江苏省境内, 上起五峰山下至鹅鼻咀, 长约91 km, 江中太平洲(扬中市)将其分为左、右两汊。2017年11月8日凌晨5：00, 该河段中扬中市三茅街道指南村附近发生窝崩险情, 至18：00时基本趋稳, 崩岸形成的窝塘呈“Ω”形状, 如图 1所示。此次窝崩共坍失主江堤440 m、岸线540 m, 坍进最大深度为190 m。本文将建立扬中河段三维水沙动力学模型, 并采用实测资料进行模型验证, 在此基础上对指南村窝崩开展模拟分析。

图  1  研究河段河势

Fig.  1  River regime of the study area

下载: 全尺寸图片 幻灯片

模型计算范围上起五峰山、下至界河口, 模拟主河道长约70 km, 包括太平洲等洲滩。基于2017年8月实测的河床地形构建数学模型, 模型采用贴体曲线非正交网格与复杂岸线贴合良好, 模型平面计算网格为1 060×240, 并在窝崩区域进行计算网格适当加密, 使网格能够较好地反映窝崩过程, 模型垂向网格则分为13层。

模型建立后, 首先采用2017年8月7—10日实测洪水大潮水文测验资料, 开展典型潮位站的潮位(水位)验证以及典型测点的垂向流速(表、底层流速)验证, 相应上游流量为43 100 m³/s左右。实测潮位测(水位)点和流速测点分布见图 1。典型测站潮位(水位)过程验证结果见图 2, 由图 2可见, 数模计算潮位过程均与原型实测潮位过程吻合较好, 包括涨落潮过程、相位和潮波变形均与原型基本一致, 各测站对应的潮位值偏差也较小, 统计结果显示最大偏差一般均在±5 cm以内, 验证结果表明河段的综合糙率为0.020~0.025。典型垂线点潮流流速过程的验证结果见图 3, 由图 3可见, 数学模型计算得到的各测流点的垂向流速过程与天然实测情况基本一致, 总体而言, 表层流速明显大于底层流速, 大部分垂线的流速值也基本与天然实测结果相符合, 相对偏差一般在10%以内。

图  2  典型测点潮位(水位)验证

Fig.  2  Comparisons between the simulated and measured water levels

下载: 全尺寸图片 幻灯片

图  3  典型测点流速验证

Fig.  3  Comparisons between the simulated and measured flow velocities

下载: 全尺寸图片 幻灯片

图 4则给出了2016年2月—2017年8月扬中河段实测冲淤分布与数学模型计算结果的对比情况, 由图 4可见, 数模计算结果与实测结果基本吻合。该时段内河道有冲有淤, 总体表现略为冲刷, 本次研究的指南村窝崩处河道主槽出现较明显的冲刷现象, 数学模型计算结果可基本反映天然河床地形的冲淤变化特征。

图  4  冲淤分布验证结果

Fig.  4  Simulated and measured results of deposition and erosion distribution

下载: 全尺寸图片 幻灯片

2.2   窝塘内水流结构特征

扬中河段为感潮河段, 指南村窝崩发生时正处于落潮过程, 水位逐渐降落。为分析指南村窝崩过程窝塘内水流结构变化特征, 以2017年11月8日窝崩发生时刻水文条件为计算边界(上游流量30 000 m³/s, 水位1.6 m), 对窝崩发展的3个阶段进行了水流结构的数值模拟分析。3个阶段模拟时除窝塘边界不同外, 其他区域地形相同, 其中第一阶段为窝崩发生初始形态, 第二阶段为窝塘已发展至一半的形态, 第三阶段是窝塘最终发展的形态。3个阶段窝塘内三维水流结构特征计算结果见图 5。

图  5  不同阶段窝塘内流速分布特征

(绿色箭头代表底层流速, 黑色箭头代表中间层流速, 红色箭头代表表层流速)

Fig.  5  Computed flow structures in the circumfluence zone during the processes of riverbank collapse

下载: 全尺寸图片 幻灯片

由图 5可见, 3个阶段的窝塘内水流结构特征相似, 均存在非常明显的回流结构。从表、中、底层流速分布来看, 窝崩初始阶段, 窝塘底层流速较表、中层流速大, 尤其是窝塘中下部这一特征更为明显(图 5(a))。对于第二阶段、第三阶段而言, 主要是窝塘尾部区域底层流速较表、中层流速大, 其夹角亦较大, 其他区域总体而言相差不大。这一水流结构特征, 为窝塘底部及岸滩坡脚淘刷提供了动力条件。

窝塘内流速最大区域基本位于窝塘内部靠里侧, 3个阶段窝塘内的最大流速具有明显的变化过程, 第一阶段(窝崩发生初始时刻), 最大流速可达0.9 m/s左右, 第二阶段(窝塘已发展至一半)时, 最大流速为0.6 m/s左右, 第三阶段(窝塘最终发展形态)时, 窝塘内最大流速为0.4 m/s左右。由此亦说明, 窝崩发生初始时刻, 窝塘内水动力较强, 窝崩发展速率较快, 当窝崩发展至后期时, 窝塘内水动力减弱, 相应的窝崩发展速率将逐渐减小, 直至达到初步平衡。

此外, 还截取第三阶段窝塘区域内上、中、下3个横断面(断面位置见图 5(c)所示), 分析断面环流分布特征, 如图 6所示。由图 6可见, 下部断面横向流速明显, 且由外部主流区流向窝塘内；中部断面窝塘内横向流速相对略小；上部断面横向流速亦相对明显, 主要由窝塘内流向外部主流区。结合图 5、图 6的水流结构可分析出:随着环流(回流)淘刷, 窝塘下部泥沙不断冲刷输移至上游, 并最终通过环流输送至外部主流区。

图  6  窝塘区域典型断面流速分布特征

Fig.  6  Computed flow structures at typical cross-sections

下载: 全尺寸图片 幻灯片

本文采用k-ε紊流模型, 模拟水流的紊动特征。第三阶段窝塘内相对紊动强度η(回流区紊动强度与外部主流区平均紊动强度之比)分布如图 7所示。由图 7可见, 整体而言, 回流区紊动强度大于外部相邻主流区, 一般为1.5~2.5倍, 平均约为2.0倍, 且窝崩口门区的紊动强度相对较大、而窝塘内部则相对较小；从表、中、底层对比来看, 最大紊动强度基本相当, 但出现的部位略有差异。

图  7  窝塘内相对紊动强度分布

Fig.  7  Computed turbulence intensity in the circumfluence zone

下载: 全尺寸图片 幻灯片

2.3   窝崩形态模拟分析与讨论

窝崩模拟时, 假设发生窝崩的上下游河岸均实施守护且不可冲, 仅考虑窝崩区域岸滩可冲刷的情况。根据现场采样分析, 岸坡下部中值粒径约0.15 mm, 上层土体的抗拉强度T₀、凝聚力c以及容重γ_b分别取值为15.1 kN/m²、12.4 kN/m²(直接快剪)和18.6 g/cm³, 厚度为2.0 m, 崩岸模拟过程中均采用此值。侧蚀冲刷系数λ采用经实测数据验证的值0.001 5。上部土体坍塌后, 一部分土体分布在床面, 另一部分泥沙则被水流带走；根据文献[18]的研究, 覆盖在床面的泥沙比例一般为30%~70%, 本文取值50%, 即假定50%的土体以源项形式进入悬沙输移, 其余50%以均匀抬高的方式平铺在坡脚网格处(1个)。模拟天然时间为13 h(即从窝崩发生时刻的凌晨5 : 00, 至基本停止发展的18 : 00), 泥沙计算时间步长为3 s。实际窝崩过程中, 窝塘内、外水体交换强烈, 水体紊动特征较强, 基于时均流速的水流挟沙力公式计算泥沙输移可能会低估河床及岸滩的冲刷过程。为分析水流挟沙力计算对最终模拟窝崩形态的影响, 本文分别采用常规挟沙力公式以及考虑回流区紊动特征影响的水沙输移方法, 对窝崩形态进行模拟。最终模拟窝崩形态结果对比如图 8所示, 图中η为回流区紊动强度与无回流区紊动强度之比, 依据相关文献[19]的试验研究成果以及前文回流区相对紊动强度分析, 本文分别取值2.0倍和3.0倍。

图  8  不同紊动强度条件下窝崩形态模拟结果对比

Fig.  8  Comparisons of predicted bank erosion shapes with different turbulence intensity

下载: 全尺寸图片 幻灯片

由图 8可见:采用常规挟沙力公式(即η=1)的窝塘形态与实测结果存在显著差异, 模拟窝塘平均崩宽约70 m, 实际平均崩宽约163 m, 侧蚀冲刷的范围明显偏小；当η=2时(即水流输沙能力增大1倍), 此时窝塘侧蚀冲刷范围较常规水流输沙公式计算范围明显增大, 模拟平均崩宽约120 m, 但仍小于实测的窝塘崩塌范围；当η=3时, 窝塘侧蚀冲刷范围明显增大, 模拟平均崩宽约160 m, 与实测窝塘平均崩宽基本相当。由此说明, 窝崩模拟时, 窝塘内水流输沙能力应充分考虑水体紊动特征带来的影响。但同时也应指出, 实测窝塘形态与模拟结果仍存在较大差异, 主要是窝塘下游侧冲刷范围模拟结果比实测结果偏大, 而上游侧冲刷范围模拟结果则比实测结果小, 一方面是因为窝塘发展形态影响因素涉及水-土等方面, 十分复杂, 尤其是实际土质条件与模拟设置的土体抗冲性不一致, 模拟过程土质设置为均匀, 而天然土质则非常复杂, 这会对模拟结果产生较大影响, 本文仅在相对简化的条件下, 对窝崩的发展过程进行了初步模拟分析；另一方面, 考虑到不同区域紊动强度的差异, 同样也会对模拟结果带来影响。此外, 窝崩的精细模拟还有待进一步深入研究, 如窝塘内回流结构是典型的非稳态过程, 其大小、结构均不断变化, 采用时均的水流模拟方法显然无法准确模拟, 可采用大涡模拟的方法模拟窝塘内的复杂水流结构；另外, 受窝塘内复杂水动力特征影响, 岸滩侵蚀机制及其泥沙输移过程均十分复杂, 还需采用更为精细的力学模式来描述和模拟窝塘内的岸滩侵蚀以及泥沙搬运过程。

3   结论

(1) 窝塘内存在非常明显的回流结构, 从表、中、底层流速分布来看, 窝崩初始阶段, 窝塘底层流速较表、中层流速大, 尤其是窝塘中下部这一特征更明显, 随着窝塘的发展, 这一特征有所减弱, 但在窝塘尾部区域底层流速仍然较表、中层流速大, 这为窝塘底部以及岸滩坡脚淘刷提供了动力条件。

(2) 窝塘内流速最大区域基本位于窝塘内部靠里侧, 窝崩发展初始时刻, 窝塘内水动力较强, 窝崩发展速率较快, 当窝崩发展至后期时, 窝塘内水动力减弱, 相应的窝崩发展速率将逐渐减小, 直至达到初步平衡。

(3) 窝塘内回流区水流挟沙力的计算对窝崩形态模拟结果存在显著影响, 采用考虑水体紊动特征影响的水流输沙公式, 可相对较好地模拟窝塘的快速冲刷过程, 因此, 窝崩模拟时, 应充分考虑窝塘内水体紊动特征带来的影响。

本文初步构建的窝崩三维数值模拟方法可较好地模拟窝塘水流结构及其冲刷过程, 为研究窝崩及其对河床演变的影响提供了一有效技术手段。