采砂河段水流泥沙数值模拟研究进展

刘亚婷，尹志刚

(长春工程学院水利与环境工程学院，长春 130012)

0 引言

当今，我国部分地区存在有法不依、执法不严的现象，非法采砂、采矿、侵占水域岸线等情形没有得到有效的治理。大量的河砂开采使得河床面下切、河床比降下降、水面线下降、边坡倒塌、河道侵蚀、渗透破坏、生态环境受损。因此，有必要应用数值计算方法模拟河道的采砂状况，从而对流域进行宏观上的掌控。不少学者通过创建物理模型和数值模型，对采砂河道演变情况进行分析。然而，给长河段大规模的河网建立物理模型时对场地要求的条件高、难度大且工作量大，因此，对采砂河段进行数值模拟分析是预测采砂影响的重要手段。

1 泥沙经验公式研究进展

20世纪初期，国外研究人员在对泥沙运移分析过程中遇到了重重困难的情况下通过现场观察和物理实验得到了基本的理论。[1]主要有DuBoy的牵引力公式(之后被Kramer修改)；爱因斯坦、Indri等开发的经典理论临界底流速度(起动速度)理论。基于这些理论，Meyer-Peter and Muller、Engelund-Hansen等开发了床沙经验公式。为了开发通用床沙公式，Garde和Albertson对其参数进行了维量分析，在这方面做出了杰出的贡献。在泥沙输移能力经验计算公式方面，Frankl、Bogardi等对其进行了研究。

20世纪末，Chang、Howard.H最早研究了采砂对河道的影响。韩其为、何明民等[2]通过研究均匀质与非均匀质泥沙运动规律及起动流速，非均匀悬移质不平衡输沙等问题，利用统计学原理，得出一系列践行于实例的泥沙运动原理，并在研究黄河泥沙问题中得到了体现；沙莫夫、列维[3]建立起以流速为主要参变量的推移质输沙率平衡公式；窦国仁、高建恩[4]以能量平衡观点建立了推移质输沙率平衡公式；黄才安、奚斌[5-6]总结出统一的水流强度指标形式(T=a1ha2Va3Ja4)，并推导出推移质输沙率公式的统一形式，如下：

未考虑起动流速时：

考虑起动流速时：

式中：K为推移质在静水中的沉降速度；Vc为起动流速；V为平均流速；d1、d2、d3、d4、d5、e1、e2、e3、e4和e5均为常数；J为能坡。为研究无黏性均匀沙的起动条件与起动标准，对国内外多种无黏性沙起动条件对应的起动概率进行了计算，分别得出起动概率值，为实际选取公式提供依据。

黄平[7]采用维里加诺夫的断面悬移质输沙率公式对广东某江修建水电站对水库河道冲淤变化进行了计算。余新明[8]采用张瑞瑾提出的水流挟沙公式与张洪武“涡团模式”得出紊动黏性系数规律，模拟计算武汉分汊河段水沙输移与冲淤演变的过程，较准确地模拟出天然分汊河道的水流泥沙特性变化。钟德钰[9]对悬移质输沙方程进行了改进，加入了反应水流曲折运动对悬移质泥沙输移影响的对流项，引入河岸变形模式。对黄河游荡型河床进行模拟研究，结果表明，方程有较好的反映岸坡变形、河势变化、洪水冲淤、散落洲滩演变的能力。郭伟健[10]运用水流方程和泥沙连续方程分析计算长江下游澄通河段水沙运动情况，得到各断面水沙变化情况，解决了大范围二维水沙计算中汊点分流分沙的计算问题。

河流数值计算中引进泥沙运动公式，弥补了数值软件中泥沙运动方程缺失的不足，加强了人们对泥沙运移特性的理解，对完善河流泥沙模拟体系起到关键性作用。为确保模拟结果的精确性与经济效益性，科学家们应进一步开发研究成本低且快捷的泥沙运动公式。

2 河流数学模拟研究进展

基于水流数值模拟的方法分析天然河道中的水力特性，主要通过使用不同方法模拟各种河型在不同水流情况下的影响来预测水流对河床演变的作用。由于河流中不仅仅存在单纯的水流运动，河床的复杂形态、水流中掺混的复杂介质都是需要研究的对象，因此，许多学者对河道的水流泥沙输移运动进行了模拟研究。

2.1 一维河流数值模拟研究

2.1.1 不含沙水流数值模拟研究

对于长河段水沙运动及河床冲淤变形计算，水流的横向运动与纵向运动相比可以近似忽略，为了简化计算，可以假定水流运动要素(流速、含沙量等)在全断面上均匀分布，建立一维水流数学模型。王领元[11]、陈德敏[12]运用MIKE11一维水动力学模型(HD模型)分别对天然河道(滞洪区与峡谷区)洪水演进及入海口河网进行了数值模拟，分析其实测水位、流量的变化情况，并进行关键断面的验证。许春东[13]、于永强[14]分别应用摩尔斯分类筛选法和LH-OAT参数敏感性分析法对河道糙率的影响因子进行局部敏感度分析，表明河道糙率在河流模拟中属于高敏感度因子，故在模拟河道水流特征时应根据河道实际情况设置糙率值。闫龙[15]应用HD模型模拟了长江上游保护区河段，分析河道内典型断面的平均河宽、深泓线、水流流速的分布规律，为保护区内鱼类栖息地生态环境特性的研究奠定了基础。付晓忠[16]应用一维水动力模型建立了吉林市城区松花江河道数值模型，分析了吉林市城市超标洪水条件下的淹没情况。

2.1.2 含沙水流数值模拟

对于一维水流泥沙数值模拟，水流及泥沙的横纵向运动比近似为0，可忽略不计，为简便计算，可假定水流泥沙运动要素在全断面上分布均匀。黄平[7]通过建立一维水沙数值模型研究在广东某江断面上修建电站和丁坝对河道的冲淤影响，研究结果表明，由于水库回水的影响，修水电站使河床绝大多数断面呈淤积状态，两江交汇处降水来沙淤积强度增大；修丁坝导致绝大多数断面呈冲刷状态，对增加航深有利。彭杨、张红武[17]建立了三峡库区即三峡大坝坝前到朱沱的非恒定一维水流泥沙数值模型，采用分流比等于分沙比的汊点分沙模式，计算得出典型断面计算水位、典型断面计算流量与实测资料符合较好，模拟出的三峡水库2003年蓄水后的排沙比与实测值接近。

目前，一维水流数学模型多用于对河道糙率的率定、水面线的计算、洪水过程的模拟与水流流速的分布情况验证。由于一维数学模型只能反映水沙要素沿水深方向上的变化，所以，对床面沿河宽变化剧烈的河湾的模拟结果可信度较低，因此，对河流的深入研究还需进行更高维度的数值模拟分析。

2.2 二维河流数值模拟研究

2.2.1 不含沙水流数值模拟研究

河道水流情况：水面宽一般远大于水深，流速等水力参数沿水深方向的变化较沿水平方向的变化要小得多，此时可将三维水沙运动的基本方程沿水深方向平均，得到沿水深方向平均二维水沙运动的基本方程。李义天[18]基于有限解析法的基本思想，用单元内插值近似地代替解析解，对长江臭盐碛河段数值模拟中的流速分布、回流范围等方面进行了模型验证。张莉[19]基于贴体正交曲线网格，应用Fluent建立了长江南通河段的二维非恒定流数学模型，分析模型对模拟天然河道水流的基本规律。

2.2.2 含沙水流数值模拟研究

余新明[8]为研究分汊河道水沙输移和冲淤演变过程，以武汉天兴洲汊道为例，采用FEM法，建立起一种适合于分汊河道的平面二维水沙数学模型。通过不同流量下流速、水位与实测值对比得出该模型能准确模拟天然分汊河道的流速水位变化、泥沙输移量、左右汊水沙比及河床冲刷淤积的演变过程，为分汊河道整治提供一定的参考价值。为研究黄河典型游荡型河道的演变进程，钟德钰、张洪武[9]建立了可模拟游荡型河流河床冲淤演变的平面二维水沙数学模型，能够较好地模拟游荡型河流中特有的散落洲滩的形成、发展、侵蚀、消亡等现象，适用于长河段、长时间的河势变化模拟。李健、杨文俊[20]模拟并讨论了在不同河段、不同方向、沙坑对河床床面变形的影响情况，通过分析沙坑对上下游流场分布情况、沙坑变形、河道代表断面面积变化和水流含沙浓度分布的规律变化得出，沙坑较下游河床的影响范围广，溯源冲刷在沙坑上游方向速度较快；沙坑的上下缘口由于较强的冲刷产生凹凸状边缘，导致局部河床加深，沙坑两侧河床发生冲刷；坑内水流中的悬移质转化为床沙，促使坑底两侧边缘产生微弱的淤积；沙坑底部发生较大的淤积，沙坑下游附近的河床冲刷幅度很快衰减。为检验垂向二维水流泥沙模型的模拟效果，邹方勇、冯卫兵[21]采用水流跨过深槽的水槽试验对水流模型进行检验，分析结果中垂向的流速、紊动黏滞系数分布、悬沙含量的垂向分布与实测资料几乎一致，表明垂向泥沙模型具有可靠性。

二维水流较一维而言更复杂，除模拟河道主流方向，还可以模拟横向水流和纵向水流的运动；主要应用于流量过程验证、水位过程验证、含沙量过程验证、河床冲淤验证等，同时体现水流对边坡岸滩以及对河床演变的影响。但在进行河床变化剧烈的河段计算的时候，二维模型不能反映出床面的变化过程，亟需进行三维河道水流的模拟。

2.3 三维水流数值模拟研究

2.3.1 不含沙水流数值模拟

由于天然水域流动得极端复杂，天然河道中的绝大多数流体运动均具有三维特性，三维数值处理主要存在水平尺度远大于垂向尺度、自由水面的非恒定运动、垂向涡黏性系数难以确定等困难。因此，对天然水域水流进行三维数值模拟是一个十分复杂而且迫切需要解决的工程问题。谭柱林[22]基于Realizablek-ε湍流模型及VOF模型模拟了辫状水槽4种工况(主槽流量为7.51 m3/s一定，主支槽宽度比B*分别为0.50、0.41、0.33与主槽流量为8.04 m3/s，B*为0.50)下的三维水流模型。分析其流场分布得出：V近水面>V近底面，辫道分流比、沙岛长宽比和水深等是分离区产生的主要因素，分离区作用下水流流速在凹岸方向有增大的趋势。吴华莉[23]应用Reynolds应力模型对正弦函数曲线衍生的连续弯道水流试验进行模拟研究，分析水位等值线、平面流速等值线、涡流强度分布得出：受弯道离心力与科氏力的影响，凸岸水位恒低于凹岸水位，各断面的水面横纵比大小不等，形成横向水面比降；弯道内底层水流流速最小，同时，弯道下游内侧区域流速减小的幅度较大，水流呈现出分离区，产生环流；连续弯道内各典型断面环流强度大小各异，水体表层和底层环流强度较大，中层的环流强度较小。

2.3.2 含沙水流数值模拟

由于沙粒自身重力作用，泥沙输运过程中泥沙浓度在垂向上存在较大差异，从而使得挟沙水流表现为密度分层的流动，加剧了流动的三维特性。

针对悬移质泥沙输移，假冬冬[24]考虑到弯道处横坡与纵坡存在对泥沙起动的影响，采用不平衡输沙方程，应用Fluent软件建立Ye(1970)的两个连续90°半径为8.53 m，中间连接处为两段长为4.27 m的直线段，进出口均有长为2.13 m的直线过渡段的弯道水槽模型和Blanckaert、Graf的单个弯道的冲淤模型，验证计算弯道环流、净冲刷、净淤积，表明除边壁处误差较大外，模拟得到的流速分布与实测基本吻合。在水流作用下，床面泥沙在湍动能的作用下被不断携带至水流中，且含沙量沿程不断扩大并逐渐形成稳定的浓度场。研究表明，两种弯道型式的水槽模拟结果均与实测值相近，呈现典型弯道“凹岸冲刷，凸岸淤积”的特点。考虑到岸滩崩塌变形对河床变形的影响，假冬冬[25]建立不同来流过程(恒定流量为0.000 63 m3/s，非恒定流其流量平均值与恒定流一致)、不同来沙情况(加粗沙：Ф=0.3 mm，加黏性细沙：Ф=0.03 mm)的三维水沙数值模型。模拟结果表明，不同流量的造床作用有很大的差别，大流量的造床效果远大于长期小流量的作用，不同流量甚至可能导致不同河型的产生；加入黏性细沙有助于散落边滩的稳固，促进弯曲河道的形成；边岸易被冲刷侵蚀，大量加入床沙，主流易散乱，有助于游荡型河道的形成与发展。模型可较好地模拟不同水沙条件下河床演变的复杂响应规律。Usman Ghani、Peter Richard Wormleaton[26]应用Fluent模拟矩形水槽中不同粒径颗粒在不同水流流速情况下的沉积效果，分析水槽断面泥沙沉积的分布情况得出：整体呈现颗粒越大沉积率越大的现象，断面主流带表现得最为明显。Robert N.Havis[27]利用HEC-6模型模拟了河道多年河床颗粒组成变化过程以及河道高程的变化过程，分别对比4种不同运移能力方程(Toffaleti、Madden、Yang、Duboys)的模拟值与实测值误差，得出最为敏感的方程Madden。Jose Oliveros-Acosta[28]利用Delft 3D对哥伦比亚的马格达莱河进行泥沙运输沉积的数值模拟，分析人为因素和天然因素导致的河床变化。在长时间水流泥沙作用下，河道呈现断面扩宽，河床深泓线加深的状况。

通过三维模拟河流水流泥沙运动，可以看出三维较一维二维水流数值模拟更形象，能直观反映出河道的横向与纵向的冲淤演变、各断面比降、各层水流流速分布等的变化情况，更有利于可视化的实现。

3 采砂坑数值模拟研究进展

天然河床由于采砂原因导致河床底部产生形状不同的沙坑，国内外不少学者对沙坑附近的流场及对河床产生的影响进行了研究。毛劲乔[29]建立各向异性的三维ASM模型(应力紊流数学模型)，对顺直河道和有边坡的矩形采沙坑附近水流形态进行模拟，计算分析水流速度矢量场及其应力分布情况得出，河床上出现采砂坑后沿主流方向产生纵向涡流，造成对上下游缘口的冲刷以及沙坑的纵向变形；坑内沿横向产生次生流，导致沙坑的横向侵蚀变形，采砂坑断面不断扩大；采砂坑底部交角处水流产生很强的剪切应力，使得采砂坑缘口部位的交角越来越小并逐渐发展成边坡，采砂坑边壁上的剪切应力不断减小，在天然河道水流运动的作用下，泥沙起动、运移、淤积，河床自行调整变形最终趋向新的平衡。为验证ASM模型的精度，毛劲乔[30]建立了形态为倒置圆锥台的采砂坑模型，模拟计算出紊流状态下沙坑的水流情况与试验变化趋势相似，这表明ASM计算精度较高，可快速预测采砂坑的影响，从而指导科学河道采砂。李健、杨文俊[31]应用VOF法设置了不同来水情况和不同沙坑几何边界条件下的模拟研究，结果表明砂坑降低了局部河道的水面线，不同情况下的砂坑影响后果不同，结果和河道与砂坑的距离，采砂坑的深度、范围以及来水来沙强度有关，影响因素效果由强到弱为砂坑的面积、沙坑的深度、来水来沙流量。李健、金中武[32]运用混合物模型对存在矩形沙坑的河道在清水和浑水两种情况下进行数值模拟，对比分析表明混合物模型可以较好地模拟河道中的水沙两相流输移现象。泥沙对流体动量的影响表现显著，且计算收敛速度明显降低，水流中的泥沙具有抑制流体内部紊动的作用。ZHANG Xiao-feng[33]应用VOF三维湍流模型，模拟矩形腔体内部不同泥沙浓度下泥沙随水流运动的分布规律。由试验结果与模拟结果得出腔体垂向泥沙浓度分布情况：层流腔体内泥沙浓度明显小于湍流；靠近上游腔体壁面泥沙浓度明显小于靠近下游腔体壁面浓度；垂向靠腔体顶部泥沙浓度最小。

前文通过在不同初始条件和边界条件下，对采砂坑进行数值模拟，将沙坑周围流场变化及冲淤变化情况对河道演变产生的作用具体的呈现出来。对顺直、弯道、分汊与游荡等典型河流不同部位进行采砂坑模拟实验，表明采砂坑的发展与河道水流流态和采砂的作用方向有关，当持续采砂的强度大，采砂量大于上游来沙补给量的时候，河床稳定将受到威胁。因此，应尽快实现科学可持续的采砂计算方法的应用。

4 数学模型的应用

目前，大多数河流泥沙运移模型具有数值性。一般而言，数值模拟得以发展，根据建模方法可以分为两组：一组基于FDM(有限差分法)，另一组基于FEM(有限元法)。这些模型也可以区分为一维、二维和三维。研究现有的数值模型发现，一维模型适用于估计或预测整个流域或河系的概况。由于计算时间方面的成本太高，三维基于FDM和FEM的模型不经常使用。二维数值模型可以为我们提供足够的局部细节和合理的成本，可以是垂直二维或水平二维。后一类也称为平均深度模型或垂直整合模型，最常用。

泥沙数学模型的发展离不开泥沙运动基本理论的研究。到目前为止，大部分泥沙输移基本理论还是建立在一维泥沙输移理论基础之上，而且现在所拥有的量测手段及对实验室和野外量测数据也基本支持现有泥沙输移理论。因此，在总量分析上，一维泥沙数学模型相对于二维模型和三维模型来说，可信度要高一些。二维和三维泥沙数学模型则在反映人们所关注的维度水沙信息上具有优势。

5 结语

5.1 结论

1)近年来，采砂问题受到广大学者的重视，河道采砂对河流本身、地表水、地下水、动植物群、航运等都有很大的影响，为控制河道采砂强度，有必要将水生态文明理念融入河道采砂的全过程，全面实行河长制，加强对河湖水域岸线的管理保护，建立健全河道采砂影响综合评估体系。

2)开发功能强大的数值模拟计算方法，借用数值方法分析采砂对河床演变的影响，建立可视化平台，实时监测演变进程，为科学可持续的开采工作打下坚实基础，实现人与河流的和谐共处。

5.2 展望

1)目前对于泥沙运动的理论在实际应用中依然存在很多的不足，亟待加强泥沙理论上的研究。

2)对于河流泥沙的数值模拟，尽量在泥沙定性的基础上实现定量。

3)由于河道各断面形态有差异，糙率也不尽相同，不同糙率值影响着水流泥沙运动与河床演变进程，应根据实际情况对河道糙率进行率定。