泥沙输移问题研究进展与展望

2019-02-14张根广王愉乐

水利与建筑工程学报 2019年4期

张根广，周双，王愉乐，刘余

(西北农林科技大学水利与建筑工程学院，陕西杨凌 712100)

床沙质、推移质及悬移质是并存于河道中的三种泥沙状态，泥沙起动、起悬是床沙进入推移运动及悬浮运动的临界状态，是探索推移质和悬移质运动规律的基础。因此，泥沙起动输移问题一直是国内外泥沙学者研究的热点问题，也是研究河床演变、水沙数值模拟、泥沙灾害防治等问题的重要基础，在实际工程中得到广泛应用。

1 泥沙起动特性

泥沙起动问题是泥沙运动基本理论中的最基础问题，也是国内外泥沙学者研究的热点问题和涉及泥沙工程中常常遇到的问题[1]。泥沙起动概念最早提出于19世纪，系统的研究则始于20世纪初，并取得了丰硕研究成果[2-3]。但由于泥沙起动问题极其复杂，所得结果并不令人满意。例如，Kramer提出的泥沙起动判别标准至今仍在沿用，致使泥沙起动标准往往因人而异，误差较大，原因在于，Kramer提出的标准是一个定性的、而非定量的标准。

床面上泥沙起动具有必然性和随机性两方面：一方面，床面上泥沙颗粒位置及作用在颗粒上的水流流速是随机的，这就造成水流强度在低于泥沙临界起动条件时，床面上仍有少量泥沙投入运动；而在水流强度高于泥沙临界起动条件时，仍有大量泥沙停留在床面。另一方面，泥沙颗粒运动又具有必然性和确定性一面。例如，当水流底速及泥沙颗粒在床面上所处位置已知时，泥沙起动又表现出确定性一面。在以往研究泥沙起动时，多采用确定性方法推求起动条件，得到的是确定性的单值关系，如Van[4]、Wright等[5]、Elhakeem等[6]。随着湍流理论的发展，学者们将湍流的随机特性应用于泥沙起动研究中，形成了随机性研究方法，如Papanicolaou等[7]。此观点认为，泥沙起动具有随机性，不存在确定的泥沙起动判别条件，由此引入泥沙起动概率问题。Einstein[8]认为泥沙起动概率P等于有利于出现冲刷外移的情况与所有可能出现情况之比，此定义包含了三种不同表述方式：(1) 作用在泥沙上的瞬时上举力大于水下重力的时间占总时间的比值；(2) 在床面任意位置，泥沙被冲刷外移的时间占总时间的比值；(3) 在任意时刻，泥沙被冲刷外移的床面面积占总床面面积的比值。但由于对近底湍流结构认识不同及测定方法的差异，学者们对近底水流特征的认识还存在较大差异，主要表现为：(1) 水流瞬时上举力服从高斯分布，如Einstein[8]、孟震等[9]；(2) 水流瞬时作用流速服从高斯分布，如韩其为[10]、Dou[11]、白玉川等[12]；(3) 水流瞬时作用流速服从对数高斯分布，如Wu等[13]、Dey等[14]；(4) 床面切应力服从对数高斯分布，如Hofland等[15]、Elhakeem等[16]，所以造成研究结果存在较大差异。此外，由于泥沙起动形式(滑动起动、滚动起动及跳跃起动)也是随机的，因此，在泥沙起动研究中，由于认识上的差异，部分学者采用了滑动起动模式(瞬时拖曳力大于摩擦阻力)作为泥沙起动判别条件，如LI等[17]；部分学者采用滚动起动模式(动力矩大于阻力矩)作为泥沙起动的判别条件，如Wu[13]、张根广等[18]；部分学者采用跳跃起动模式(瞬时上举力大于水下重力)作为判别条件，如Einstein[8]、Cheng等[19]、Dey等[14]。但由于泥沙运动极其复杂，在实际观察中，通常很难准确地将三种起动模式区分开来，一些学者只好采用概化理论公式推导三种模式的临界水流强度，如Dey等[14]。

床面泥沙三维分布特性极其复杂，影响因素众多，概括起来可表述为：表层泥沙密实性、泥沙非均匀性、绝对暴露度及相对暴露度等参数。早在1950年，Einstein[8]就发现了采用均匀沙代替非均匀沙引起的推移质输沙率的误差，并通过引入相对暴露度系数试图修正非均匀沙所引起的上举力偏差。之后，Paintal[20]通过分析理论计算床面和平均床面高程，首次定义了暴露度新概念。韩其为[21]在对泥沙起动流速研究中，基于床面泥沙颗粒相对位置，提出了概念清晰且便于应用的绝对暴露度和相对暴露度概念。刘兴年等[22]通过对非均匀沙的暴露高度进行实验测量，给出了暴露高度与粒径间的关系。近些年，学者们分别从理论分析及实验观测两个方面，对床面泥沙颗粒位置进行了深入研究，何文社等[23]提出了暴露度定义的新方法，并给出了对应的计算公式，得到了基于等效粒径的非均匀沙暴露度系数取值范围。杨奉广等[24]采用Cheng等[25]的床面颗粒间相对步长概率密度分布的研究成果，结合床面泥沙暴露度与步长之间的几何关系，得到了河床泥沙颗粒相对暴露度的概率密度分布函数；孟震等[26-27]对二维、三维泥沙颗粒相对隐蔽度进行了初步分析，并对其物理意义进行了讨论；白玉川等[28]提出了双向暴露度的新概念，重点研究了垂向暴露度和纵向暴露度对泥沙起动流速的影响；周双等[29]、邢茹等[30]通过自行设计的实验装置，测量分析了随机堆放条件下均匀玻璃球及河卵石纵向水平间距、暴露角、绝对暴露度及相对暴露度，得到了玻璃球及河卵石相对暴露度的概率密度均服从偏正态分布。需要指出的是，以上研究虽然一定程度上揭示了泥沙位置的分布特性，但大部分学者多考虑的是均匀沙，这与实际河流中泥沙为非均匀沙、且长期受水流作用的影响还有一定出入。

总的说来，对于泥沙起动的研究，还未脱离确定性方法的老路。在将来的研究中，还应多考虑采用统计分析方法来描述泥沙起动状态，既要考虑水流的随机性，也要考虑床面泥沙非均匀性、随机分布、泥沙颗粒形状的不规则性及颗粒大小之间的隐暴效应等因素。

2 推移质运动规律研究

早在19世纪末期，法国Duboys基于水流剪应力，开创性地提出了推移质泥沙输移理论，之后众多学者对这一问题进行了深入研究[31-37]。从研究方法上来说，推移质运动的研究方法可划分为实验研究方法与理论研究方法。实验研究方法以Gilbert[31]及Meyer-Peter等[32]为代表，理论研究方法的途径则较多，且各有千秋。Gilbert在1909年—1914年，采用了9种沙样进行了892组试验，开创了水槽试验的先河，标志着现代泥沙运动力学的兴起，其研究成果至今仍在被引用；他们在孤立单因素分析基础上，通过系统的推移质输沙率与单因素、多因素的试验，逐步建立和完善了以水流切应力为水流强度指标的推移质输沙率公式。在理论研究方面，有基于流速为主要参变量建立的推移质输沙率公式，如沙莫夫公式、冈恰洛夫公式，列维公式等；有基于拖曳力为主要参数建立的推移质输沙率公式，如耶格阿扎罗夫公式，恩格隆公式，阿克斯-怀特公式等；有基于能量平衡观点建立的推移质输沙率公式，如拜格诺公式，杨志达公式等；有基于统计方法建立的推移质输沙率公式，如爱因斯坦公式，韩其为公式等；有基于沙波运动建立的推移质输沙率公式，如shionhara公式等[33]；近几年来，张根广课题组[17-18,29-30,34-39]将水流随机性及泥沙相对位置随机性同时引入泥沙起动输移中，从理论上推导出全概率泥沙起动及输沙率公式。随着非线性科学的发展，一些研究者试图从非线性角度对推移质输沙率进行研究，例如，白玉川等[40]通过分析推移质运动过程的非线性动力学特性，证明了床面形态的演化决定于推移质运动过程的突变特性；杨具瑞等[41]基于尖点突变理论，得到了推移质输沙率公式；在推移质泥沙运动的分形和自组织规律研究方面，姚令侃等[42]基于非均匀沙水槽实验，建立了推移质输沙率与床面泥沙颗粒间的双曲关系，但这种分布的动力学机制尚需要进一步研究；汪富泉[43]通过理论分析和实验数据论证，论述了推移质运动的自组织临界性和分形列维稳定分布的物理机制。综上分析，上述研究成果大多数均建立在均匀沙的基础上，即使部分成果是基于非均匀沙研究，但由于分析中均采用的是泥沙中值粒径或平均粒径，因此不算真正意义上的非均匀沙输沙率研究，而完全基于非均匀沙考虑的推移质输沙率计算方法有四类[44]：直接分组计算法、修正剪切力法、床沙分组计算法、输沙率级配法。直接分组计算法采用床面泥沙分组计算输沙率方法，得到不同粒径组的泥沙输移率，最终累计求和得到总输沙率，如刘兴年等[45]、孙志林等[46]、陈有华等[47]、Cao等[48]。修正剪切力法是通过引入剪切应力修正系数，得到床面不同粒径组的床面剪切力，将修正后的床面剪切应力引入到均匀推移质输沙率中，得到不同粒径组的推移质输沙率，如Parker等[49]、Wilcock等[50]、Rickenmann等[51]、孙东坡等[52]。床沙分组计算法，首先要确定床面总的可能输沙率，其次要确定床面泥沙级配，而后将床面总的输沙率与不同粒径组在床沙中的权重乘积作为床面分组输沙率。该方法的计算精度受粒径分组数目的影响较大，现有的计算公式精度均还有待提高。输沙率级配法类似于床沙分组计算法，首先要确定床面总的可能输沙率，其次是确定推移质运动泥沙级配，而后将床面总的可能输沙率根据推移质运动泥沙级配进行分配，从而确定各粒径组的输沙率。该方法有效的避免了床沙分组计算法中分组数目对总输沙率计算结果的影响，提高了公式的预测精度，如窦国仁等[53]、Hsu等[54]、Tayfur等[55]、Bombar等[56]。

综上可见，无论是均匀沙还是非均匀沙，目前的研究成果多以单颗粒泥沙及单一运动状态为研究对象。而在实际河流中，随着水流强度的增强，床面泥沙并不是以单颗粒泥沙运动和单一运动状态运动，而是以多种运动形式(滑动、滚动、跳跃等)和群体运动状态输移。惠遇甲等[57]、Hu等[58]采用高速摄像技术对床面附近颗粒的运动状态进行了观测研究，初步得到泥沙滑动、滚动、跃移以及悬移等各种运动形式所占的比例，开创了研究群体泥沙运动状态及多泥沙运动形式的先河；Gao[59]研究了泥沙滚动和跳跃之间的转化、单层和多层运动的转换。

在泥沙运动特性研究中，早期主要以泥沙受力结构为基础，通过构建显式力学平衡方程，求解得到推移质运动速度与水流强度之间的关系。随着近代计算技术的发展，使得求解复杂的微分方程成为可能，因此，基于经典力学，结合近代流体力学、概率论、随机过程等学科，构建基于单颗粒受力模型、考虑颗粒位置随机性或受力随机性的微分方程，通过数值求解泥沙颗粒运动速度取得一定进展，如Lee等[60]、白玉川等[61]、徐海珏等[62]、范念念[63]。

在实验研究方面，学者一方面结合最新的理论进展，尝试采用多影响因素进行实验研究；另一方面不断提高测试手段，例如高速摄影技术、图像识别技术、粒子示踪等。惠遇甲[53]利用光学补偿照相机，刻画了单个泥沙颗粒的跃移轨迹，对颗粒的速度、加速度、旋转、受力等问题开展了研究；唐立模等[64]基于三维粒子示踪测速技术，观测了10种水沙条件下推移质颗粒的三维运动情况，统计分析了推移质颗粒的三维平均运动规律；余国安等[65]通过专门设计制作的推移质运动视频观测采集系统和河床结构发育程度测量装置，对几种经典床面条件的山区河流推移质颗粒强度运动状态进行了观测分析；Julien等[66]通过对光滑床面上的不同大小、不同形状和密度的泥沙颗粒的平均推移运动速度进行了观测，拟合得到了推移质运动速度计算公式；胡浩等[67]采用ADCP定点观测方法，通过对长江河口南港和北港粉砂、细砂河床的推移质运动速度进行研究后指出，只有在床沙发生普遍运动时，ADCP观测法才能进行有效的测量；Amir等[68]采用摄像粒子图像测速系统及脉动压力传感系统，测量了床面颗粒跃移参数及受力变化；等等。与均匀沙相比，非均匀推移质运动更加复杂，主要表现为不同床面形态下，河床粗糙度对颗粒运动的间歇性和变化程度产生直接影响。Nicholas[69]指出，颗粒顶部湍流强度将会随着颗粒糙率的增加而增大，这一变化将引起颗粒速度进一步加大；Ramesh等[70]研究了水力过渡区内泥沙运动速度，得到了泥沙运动速度的分布规律；许琳娟[71]在水槽中开展了非均匀推移质运动特性的实验，利用数字图像处理技术，获得泥沙位置及运动轨迹，对泥沙运动速度、走停时间等参数进行了研究。

总之，在推移质运动及输移方面虽然取得了丰富的研究成果，但大多成果主要基于单颗粒泥沙运动特性，很少有从群体泥沙颗粒运动角度去探索推移质输移规律，尤其是从群体颗粒内部颗粒间的作用力(例如颗粒碰撞力)角度去考虑。就是对于推移质运动特性的认识，仍停留在对相关现象的定性描述，缺乏对这种现象产生原因的认识。究其原因是多方面的，既有学科面窄，无法与基础科学、尖端技术等相比的原因，也有学科本身问题复杂、难度大，而测量仪器不能满足科研需求的原因。

3 悬移质运动规律研究

泥沙浓度分布是河流动力学中基本理论研究的重要内容之一。自1930年代，Rouse[72]基于紊动扩散理论，首次建立了泥沙浓度分布公式以来，至今己有多种不同的理论模型，如混合理论、能量理论、相似理论及随机理论。而在众多研究成果中，又以紊动扩散理论居多，其中扩散系数的选取是导致浓度分布公式繁多的主要原因。倪晋仁等[73]基于垂向脉动速度概率分布特性，推导出了悬移质浓度垂线分布的统一公式及悬沙紊动扩散系数表达式，并从理论上揭示了这些理论的出发点虽然不同，但最终都能得到与扩散方程类似结构形式，只是扩散系数略有不同；并进一步分析提出，含沙量沿垂线分布包括两种类型：(1) 泥沙浓度沿垂线分布规律为从水面到河底由小逐渐增大，在近壁区达到最大值，然后再由最大值迅速减小为零，称为Ⅰ型(普通型)浓度分布；(2) 自水面到河底单调增加称为Ⅱ型[74]浓度分布。然而，采用传统的泥沙扩散理论很难解释室内试验和原型观测到的Ⅰ型泥沙浓度分布，因此，学者门试图寻求其它理论和方法解释这一现象。如Wang等[75]基于固液两相流动力学理论分析后指出，Ⅰ型泥沙浓度分布规律的形成原因，主要来自于升力的影响。刘大有[76-77]早期基于两相流的紊动扩散模型研究发现，传统的紊动扩散模型不能很好的描述颗粒脉动强度梯度引发的扩散；尤其是在固相颗粒脉动强度变化较大的区域，基于菲克第一扩散定律的泥沙扩散理论更是存在明显的缺陷，并定性解释了Ⅰ型固体浓度分布的成因；之后基于泥沙悬浮运动机理认为，悬浮泥沙颗粒在河流或水平管道中，除受到重力和浮力作用外，还将可能受到流体中其他升力和泥沙颗粒脉动强度梯度的作用。傅旭东等[78]基于流体动力学理论，修正了传统的泥沙扩散方程，并对影响含沙水流泥沙浓度分布类型、形成机理及扩散系数的影响因素进行了定量分析。倪志辉等[79]基于非线性方法中分形标度，从流体内部结构分析入手，探讨了挟沙水流的含沙量分布。

1950年代，Einstein等[80]基于床沙质、推移质和悬移质三者之间的内在联系，构建了水流挟沙力理论公式。Bagnold[81]基于水流功率理论，结合推移质、悬移质的研究结果，推导出了一个包含推移质和悬移质的全沙输沙率公式。可见，上述两个公式概念明确，具有较强的理论基础，但其适用性与经验性公式相比，还有一定的差距。例如以悬移质泥沙输移为主体的平原河流，推移质泥沙运动可以忽略的情况下，经验性公式更具适用性，类似的公式有武汉水利电力学院学者[82]公式、Velikanov[83]公式及范家骅[84]公式。综上可见，上述几个学者的研究都是基于推移质和悬移质具有不同运动形式的基础上，在具体研究时区别对待，沙玉清[85]则认为，推移质运动和悬移质运动都是泥沙一种运动形式，在本质上没有什么不同，并基于这一思想，构建了一个更具普遍性及适用性的公式。

高含沙水流具有与普通挟沙水流不同的运动特点和冲淤规律。自1980年代以来，高含沙水流的挟沙能力问题，得到了国内外学者的普遍关注，并确定丰富的研究成果[86]。曹如轩[87]结合陕西省水利科学研究所室内试验资料，通过对高含沙水流的浑水粘滞性、泥沙沉降特性、高含沙水流极限粒径的概念及输沙特性等因素分析研究，建立了高含沙水流挟沙力公式；张红武[88]基于泥沙悬浮能量来自于水流运动动能、水流能量消耗与泥沙悬浮功之间必存在着内在联系的考虑，通过分析研究泥沙含量对卡门常数与颗粒沉速的影响，以及浑水粘滞系、水力泥沙因子对高含沙水流挟沙能力的影响，提出了一个半经验、半理论的高含沙水流挟沙力公式；吴保生等[89]采用回归分析方法，分析研究了水流挟沙力的内在机理及影响因素，构建了基于浑水黏性系数、泥沙沉速的高含沙水流挟沙力经验公式；费祥俊[90]基于黄河泥沙悬浮指标推导出黄河干流不冲不淤流速，并在建立黄河干流挟沙能力计算中，计入黄河床沙级配及河道断面形态特征的影响，得到了一个适合黄河挟沙能力计算的经验性公式；舒安平等[91]基于泥沙悬浮和水流紊动之间的联系，借鉴两相流的相关研究成果和大量实测资料，构建了高含沙水流挟沙能力公式。总之，目前在悬移质浓度分布、水流挟沙力、高含沙水流等方面取得了不少研究成果，但总的来说，水沙两相流的理论发展还不够成熟，现有的理论研究成果或公式还存在较大缺陷。因此，在将来很长时间内，还要重点研究低浓度泥沙颗粒与水流间的相互作用，揭示颗粒悬浮、扩散与交换机理；在此基础上，研究高含沙对水流运动特性、能量耗散特性的影响，探求高含沙水流与一般含沙水流在水沙输移方面的区别与联系。这就要求研究者不断加强基本理论的研究与探索，不断提高试验观测技术水平，着力实现跨学科、跨领域合作，最终建立一个适合含高含沙水流的悬移质浓度统一分布形式与水流挟沙力公式，在体系上完善悬移质输移机理。

中国北方灌区多属于干旱地区，灌溉水源多为渠灌引水。尽管灌渠挟沙水流的输移规律的研究可借鉴明渠挟沙水流的运动理论，但仍与天然河道存在一定差异，一方面，灌溉渠道边界条件规整，便于研究与应用；另一方面，由于不同灌区(或同一灌区不同时段)来沙特性、水力条件、以及管理方式的差异，灌渠泥沙的输移规律也表现出特殊性。相较天然河道而言，灌区渠道泥沙淤积给灌区正常运行带来了严重问题。因此，研究者对灌区泥沙运动特性及泥沙淤积成因(泥沙淤积与来水来沙间的机理)进行了研究。早在1987年，人民胜利渠灌区运用之初，河南人民胜利渠泥沙专题研究组[92]就灌区水沙分配规律和泥沙淤积分布规律开展了细致的研究；随着工程的不断运行、改造，以及黄河水沙条件的变化，相关研究不断深入，王延贵等[93-94]以渠道水流含沙量垂线分布和输沙能力为基础，对引水分沙的特性及其对渠道冲淤的影响进行了初步分析，之后，在大量水沙实测资料基础上，对泥沙起动规律、含沙量分布规律及挟沙能力进行深入分析，指出渠道挟沙能力主要受渠道糙率、渠道比降、横断面形状、来沙级配和引水流量的影响。此外，李春涛等[95]、史红玲等[96]对位山引黄灌区的泥沙淤积原因也进行了相关研究，其结论与王延贵等[94]的结论基本一致。黄河水沙的总体特征是粒径小，含沙量高，而国内部分渠道所引水沙与黄河水沙并不相同，故输移规律也并不一致。山溪性河流的引水问题与黄河的引水问题的最大区别在于，除了悬移质外，对推移质的输移规律同样需要给予关注[97]。

为了兼顾灌区浑水资源的可持续利用及防淤防堵两大目标，部分研究者以泾惠渠、渭惠渠灌区为例，对渠系节点处的泥沙淤积特征以及对引水分流特性开展研究。王延召等[98]对渠系节点区域悬移质淤积与支渠引水分流间的关系进行了研究，指出回流区水流挟沙能力减小，泥沙在口门前形成拦门沙坎，是产生淤积的主要因素。同时，为了实现灌区水沙资源的优化配置和泥沙运用，研究者们进行了大量的研究，史红玲等[99]构建了基于入渠水沙资源分配目标的灌区配置能力指标的表达式，并对指标值进行了量化；胡健等[100]指出灌区泥沙合理配置与利用的关键是非均匀沙的输送机理，进而对非均匀沙不同粒径组泥沙的起动与止动、沉降与悬浮的差异进行了研究，得到分组沙的输沙能力、沿程衰减规律和上限平衡含沙量的统计特征。