冲积河流糙率由来与计算方法研究

张红武，张罗号，彭 昊，蔡蓉蓉

（1.清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084；

2.河海大学 水利水电学院，江苏 南京 210098；3.中国科学院 地理科学与资源研究所，北京 100101）

1 明渠均匀流流速公式研究回顾

在“动能”和“势能”术语还无定义之前的1738年，瑞士数学家丹尼尔·伯努利（Daniel Bernoulli，1700—1782）在出版的专著[1]中描述了流体的压力与流速的关系，阐述了“活动力”的守恒，尽管其研究结果跟现在众所周知的伯努利方程有相当距离，但仍为伦哈德·尤勒（Leonhard Euler，1707～1783）等建立水流能量转化与守恒定律提供了条件[2]，因此后来的学者把水流能量方程以“伯努利”冠名[2]，并不完全是个误解。无论如何，水流能量方程的研究进展无疑为以后建立明渠流速的计算方法打下了基础。多年后的1768年，为从伊沃特河修建一条引水渠到巴黎市区来解决供水紧缺状况，市政主管部门委托当地水工建筑学校首任董事长佩罗内特（Perronet）与教师安托万·谢才（Antoine Chezy，1718—1798）开展引水工程设计。谢才分工负责设计明渠的断面尺寸[2]，为找到计算方法，他按照均匀流水流重力与阻力应该相等的物理条件，假定在一定时间里运动粒子所受阻力（亦即明渠壁面的切应力）与断面平均流速的平方成正比，同样也与湿周长度成正比，1769年9月中旬前推导出半理论的关系式[2-3]。教科书常认为谢才公式是此时提出的[4-5]，其实1769年建立的仅是公式的初步形式，至于公式是1755年[6]或1768年[7]提出的也不正确。另外，教科书[5，7-10]无视公式建立在正确物理图形上的事实，称谢才公式为经验公式颇不妥当。当年为检验公式的系数[2-3]，谢才在风和日丽的9月23日与10月7日，分别选择科尔帕里特（Courpalet）运河和塞纳河（Senie）的顺直段进行了两组近似恒定均匀流河渠观测试验，通过蜡球浮标测定水面流速[2-3]。有教科书称“1769年谢才对河渠均匀流进行试验研究后”[5，8～10]才建立公式，也混淆了他的研究路线。谢才利用试验资料对初步公式反算流速系数发现，两组试验的结果差距较大（第一组约为第二组的71%）[2]，自认为是流速测验误差过大所致，不过他仍通过对第一组观测结果修正，推断出拟建引水渠过流能力“绰绰有余”。按谢才当年学生巴伦·里奇·普罗尼（Baron Riche de Prony，1755～1839，法国杰出水利科学家）1804年备忘录的记载，谢才1775年就确立了他的公式[2-3]，在1776年的备忘录上才将自己的流速公式简化成如下形式[2-3]：

式中：V为平均流速；J为能坡或比降；R为水力半径；C 在1897年上式公布后被称之为“谢才系数”，并在曼宁（Robert Manning，1816—1897）公式出现并经同上式比较后，被认为是河床糙率n及水力半径R的函数。

为试图将系数C 确定为常数，采用上述两组试验修正数据反求式（1）中的C 大致为57.3[2]（英制单位，国际单位或公制单位相当于31.6 m0.5/s），故有些教科书讲“谢才曾认为系数C是个常数，并取50 m0.5/s”[7～10]实际是个误解，何况直到19世纪上半叶各国流行的几个与式（1）形式相同、C为50 m0.5/s 左右常数的公式，都是以其他学者名字命名的[2]。尽管供水系统设计负责人佩罗内特在伊沃特引水工程设计中采用了谢才的试验结果，但在设计报告书中并没有包括他有关式（1）的分析内容。谢才在任教的学校退休后生活穷困潦倒，待往日学生普罗尼帮他1797年担任本校董事长的第二年即辞世[2]，也同他没认识到自己成果很有价值有关，使式（1）长期没有公布。多年后的1803年与1804年，其学生吉拉德（Girard，1765—1836）和普罗尼先后在备忘录中提及搁在档案室的谢才报告[2]，但没引起法国同行的重视。直到1897年，德国土木工程师克莱门斯·赫谢尔（Clemens Herschel）将这份“积满灰尘”的报告，整理并翻译成英文在美国发表[3]。此时流行一个多世纪的是同样采用“水流阻力同流速的平方成正比”假设的流速公式，在1786年是法国军工工程师杜博（Pierre Louis Georges Du Buat，1738—1809）发表的[11]。谢才与杜博所处年代“人人都谈论水力学”[12]，“使水流运动的力应等于由于流体黏性引起的阻力和接触面的摩擦力之和”这一概念被人们普遍接受，Du Buat公式结构繁杂，但他是首次发表的用物理分析方法描述其原理的学者。

巴辛（Hanri Emile Bazin，1829—1917）1854年成为达西（Henry Philibert Grspard Darcy，1803—1858）的助手后，共同开展了大量室内实验，依此1865年提出了Darcy和Bazin公式[13]：

将上式与式（1）比较，即可得知（RJ）0.5之前有一个类似的系数C（当时谢才不为人所知，故称C为“流速系数”），引起了学术界关注，但式（2）中两个待定系数限制了实用价值。巴辛认识到该问题后，引入粗糙度系数n1（比常见的糙率n大很多）1897年给出了C的新公式（与谢才公式同年被发表，后来我国学者习称C为谢才系数）[14]（英制单位）：

1869年瑞士学者Ganguillet及Kutter 首次引入糙率n（又称粗糙系数[15]，有些教科书[4]受现状英文习称“Manning’s roughness coefficient”影响，将其称为“曼宁糙率系数”或简称为“曼宁系数”，实不妥当），给出计算C的如下公式（C的单位m0.5/s）[15]：

上式跟式（3）都存在量纲不和谐的问题，尤其在R＞1、＜1、等于1 m 三种情况下，C随J的变化关系相应呈正比、反比、无关的三种规律，其表现未必合理[16]，特别是跟冲积河流冲淤变形同J与摩阻特性变化趋势不符。尽管如此，式（4）问世后很快在欧美水利工程界广泛接受[17]，大量积累出不同边界材料的糙率n取值表[16]，由此“通常给出令人满意的结果”[18]。其实，已是水利专家和技术官员的曼宁，也是该公式的应用者，他认为这是个“极有价值的公式，普遍受到青睐”[2]。于是，年愈古稀的曼宁以此为基础建立了形式简单的流速经验公式，1890年发表论文的公式形式为[19]：

式中：c为流动阻力系数，文献[17]认为是“经后人修改为”以下常见形式（即c=1/n）：

从曼宁论文只有式（5）而没有式（6）的情况看，文献[17]的推断是有道理的。不过，从曼宁在研究式（4）时早已深知糙率n对计算流速的重要作用等情况判断，曼宁发表式（5）前已确定出c 同n 呈反比的关系，否则他就无法根据实测资料确定出式（5）中R的指数为2/3（J的指数为1/2 当时已成共识）。

由于式（5）的因次不和谐，且像R的指数2/3 在当时计算是相当累赘的，故曼宁对式（5）并不感兴趣，也许跟此前几十年内已出现不少和式（5）形式相同的公式有关。例如，菲利普·加尔帕·高克勒（Phlippe Gaspard Gauckler）1867年提出的两个公式[20]中的第二个，即J＞0.0007的公式与式（5）形式相同；再如，1881年Hagen 也得出相同的公式[21]。以至于曼宁在1889年12月4日交给爱尔兰土木工程院的论文之前言中坦承[19，21]：“尽管这个公式是作者在1885年独立得出的，但坎宁安（Allan Cunning⁃ham，R.E.）在他的论文《当代水力学实验》中已经提到过；近来Hagen 博士根据库特的实验结果，用最小二乘方法也导出了公式V=cR2/3J1/2”。

于是，曼宁1885～1894年继续通过考虑大气压力的影响，试图找到一个既能用于测算管道流量又能用于计算明渠流速的通用公式[2，21]。不过，这个新公式提出之前，1891年即推荐了式（5）[22]，几年后同行纷纷称此式为“曼宁公式”，且1911年Buckley还在工程中进行了应用[21]。在Ganguillet及Kutter公式积累的丰富糙率n数值表基础上，工程界经过对曼宁公式的应用进一步充实了n取值图表；且发现只要n值选择恰当，曼宁公式可满足一般设计要求[18，23]。正是这简单而重要的理由，使这个缺乏理论基础、量纲不和谐的式（6），成为现在水力学中运用最为广泛的曼宁公式[23]。至于他花费9年时间提出的通用公式，因没涉及同糙率类似的系数而无法在千变万化的实际中运用。到了19世纪末，包括式（4）在内的很多水流阻力公式逐渐不被使用。但也不像有些文献所讲“均遭淘汰”，甚至称是“去伪存真，去粗取精”。其理由一是后人对当年成果挖掘不够（例如，现不常见的Ganguillet及Kutter公式，在某种条件下就可能减少采用曼宁公式出现的失误）；二是现有方法并不完美[24-25]，出现失误后又都归结到糙率取值不当问题[7-8，10]。

从现有教科书看到，几乎都称下式为“曼宁公式”[4-5，7-10，17]：

其实上式只是后人联解式（1）与式（5）得出计算谢才系数C的表达式，毕竟曼宁辞世那年谢才公式被发表后才有“谢才系数”之称，故不应称式（7）为“曼宁公式”，而应称为“基于曼宁公式的谢才系数计算式”。

19世纪上半叶学术交流不够，很少人注意到文献[3]，各国流行一些与式（1）形式相同、系数C为50 m0.5/s 左右常数的公式（例如德国Eytelwein公式与意大利Tadini公式等）。若没有式（4）与式（7）或者式（3）之类的“流速系数”计算式，只有等到20世纪30年代尼库拉兹[26]、蔡克士大等学者的试验丰富达西1858年的管道阻力资料[27]，使H.Darcy-J.Weisbach公式[14]中的沿程阻力系数λ有相应取值后，由确定谢才系数[17]，式（1）方有使用价值。正是人们逐渐掌握了谢才系数C同n、λ之间相互换算关系，才使式（1）成为水力学十分著名的谢才公式。附带指出，水力学界学者知道，因为沿程阻力系数λ是无量纲的，H.Darcy-J.Weisbach 平均流速公式比曼宁公式严谨，但应用远没有后者广泛，其原因是λ 同水深、边壁粗糙程度的关系过于复杂。

在前苏联，1925年巴普洛夫斯基分析所积累的丰富试验资料，指出式（6）水力半径R的指数取为常数，“往往造成重大的错误”，为此建立了谢才系数C 同糙率n的如下关系式[17]：

式中指数y和水力半径R、糙率n有一定关系，即：

计算表明，对于糙率小于0.03的沙质河床而言，上式与式（6）差异不大[17]。

人们称式（6）为“谢才公式”也是不严谨的，至多可称为“谢才-曼宁公式”。实际上，近一个世纪河道与渠槽的水力计算，采用的是曼宁公式或所谓“谢才-曼宁公式”，而基本没有应用取C为常数的谢才公式（1）的实例。由此说明，糙率n对于水力学及河流研究十分重要。

2 动床糙率计算方法研究现状简述

由上述可知，1869年在甘固利特及库特公式中首次出现的糙率n，并在广泛应用过程中积累出糙率n数值表，使结构简便的曼宁公式一经推荐即脱颖而出。然而，对于糙率n的定义，人们并不那么关注。一般可称其为“反映壁面粗糙程度的系数”[5]。由于曼宁公式是量纲不和谐的经验公式，对于糙率n一般不写单位或量纲，既不需要像有的教科书[8，10]那样强调糙率为“无量纲系数”，也不能像有的教科书[4]称糙率“具有量纲L1/2T”，以免在国际单位同英制单位等相互换算时出现麻烦，尤其可能使长期工程经验积累的糙率n数据表变乱。

对于天然河道来讲，糙率n定义就显得更为复杂而重要[24-25]。我们以其涉及的主要方面，在不考虑水草等情况下，试将河流糙率n定义为“表征河床湿周形状规则性、粗糙程度、河底形态与水沙情况对水流阻滞作用的综合性影响系数”。在实际的水力计算中，由于目前糙率图表给出的条件是定性的，对应的n 又多是范围值，因此糙率n取值因人为性较大而难以确切。设计人员习惯根据沿程水位实测点据或设计需要，分若干河段反调糙率，戏称“糙率是个筐，什么都往里面装”，调出的各个河段糙率值忽大忽小，相当于基本概念表现往往超常。其实，这主要是没考虑流量变化[6]、床面泥沙运动[28]、含沙量存在[29]、河道平面与断面形态[25]等方面对水流摩阻特性的影响所致。

水力学领域对尼库拉兹1933年利用管道系统开展的均匀沙粒阻力试验[26]，颇为熟悉。其实，在此之前不少学者已将糙率直接同泥沙粒径D 相联系，因受上述式（7）这个基于曼宁公式的谢才系数计算形式的影响，自动想到粒径D的指数为1/6，亦即沙粒糙率与D1/6成正比，因此，1923年Strickler.A根据实测资料提出如下沙粒糙率关系式[30]：

式中D 常取用床沙组成中以重量百分比计65%较之为小的粒径D65，mm。

对于粗沙平整河床，上式基本符合实际。但对于细沙河床，计算值较实际糙率偏小不少[24]。为此，张有龄1939年通过水槽试验，将上式中的系数0.015 调整为0.0166[31]（对于细沙河床，一般文献认为D 取用的是中值粒径D50，与D65的近似关系[32]为D50≈0.82D65）。后人将粒径D以mm 计并取其式中系数的倒数，符号记为A 但却只称为“系数”，由于该系数A= D1/6/n，始称之为摩阻系数[32]。鉴于该系数已被不少学者作为参变量加以研究[33-34]，再参考英文有时称为“empirical roughness parameter”，建议称之为“摩阻参数”。则将沙粒糙率的计算公式表示为：

上述张有龄公式，A相当于19；Strickler公式A相当于21，但他对于莱茵河卵石河床，取A=24。将上式代入式（6），得

对于静平床而言，A的取值与床沙粒配、形状及排列状况有关，故不同学者根据具体的床面情况，调整了摩阻参数及代表粒径取值。如果代表粒径取D65，摩阻参数符号也可记为A65，一般可采用下式计算[32]：

式中D0为避免量纲不和谐而引入的参考粒径，取D0=1m。

河床出现冲淤时床沙粒径随之粗化和细化，糙率相应增减，看出式（10）反映沙粒的肤面阻力变化是符合实际的，因此，这类关系式也属于动床糙率公式。钱宁1955年归国后同早期归国的黄河专家麦乔威团队积极合作，在黄河泥沙研究领域取得了大量成果。其中他们以式（10）为基础[33]，选取黄河下游花园口等水文站资料，分析了摩阻参数A（代表粒径取D65）与Einstein 研究推移质运动时提出的水力参变数ψ（简单换算后基本与Shields 数类同）的关系，使通过摩阻参数A计算糙率的条件不仅包括了水流与床沙粒径条件，而且也不自觉地反映了含沙量的影响（尽管没有将含沙量作为影响因子，但选用的都是浑水资料），条件远超过了沙粒糙率限定的静平床。鉴于他们给出接近动平床时的糙率公式与张有龄公式（A=19）相同，可以认为在研究黄河问题时，不能完全视张有龄公式为清水静平床沙粒糙率公式，否则就不能解释由此式计算的黄河下游糙率为何比玻璃糙率还小不少的现象。此外，马睿近两年以此为基础开展的研究[35]，定量描述了更高强度水流条件下摩阻参数变化规律，使多沙河流输沙能力研究同河床阻力变化有机联系起来。

这方面更系统研究的是1963年李昌华及刘建民的成果，他们分析黄河、长江，赣江等天然河渠的资料，点绘出摩阻参数A（代表粒径取床沙中值粒径D50）与相对流速V/Vc之间的关系图（V为平均流速；Vc为泥沙起动流速），对于长江，上述沙粒糙率公式中D50的指数为1/6，对于黄河、赣江，D50的指数取1/5，由此可包含由静平床到不同河床形态的相关过程。再将对数流速分布公式中的ks表示为D50与无量纲系数α的乘积形式，以α与相对流速V/Vc的函数关系，来反映水流强度引起床面形态变化对水流阻力影响[34]。

Engelund[36]等西方学者习以对数流速公式为基础，通过理论假定与推导，给出不同的无量纲参数及其函数形式，再整理实测资料（多是水槽试验资料）得到经验曲线，即可在已知R、J及床沙粒配条件下，试算求出V。这类方法过程繁琐，且根据实测资料给出的曲线可靠性与代表性都差，即使一些学者历经半个多世纪的修正，但使用价值依然有限。

沿程阻力系数λ与糙率n均是无法直接测量出来的物理量，前者是由λ=8 gRJ/(V2)的关系，后者则用曼宁公式（6），利用测验的V、J、R 反求的。由于水文站对J采用的是测验误差大的水面比降值，故反求的糙率λ、n的精度很低。对于宽滩窄槽的河流（H≈R），考虑漫滩后的断面平均流速V会变小不少，平均水深H 变小更多，采用主河槽测量的水面比降作为J，再用公式反求的λ、n值精度更低，还常出现比自然界最小糙率还小很多的异常现象[24]。例如，利用我们开展的黄河封丘河段动床模型试验结果可发现问题。该模型按照黄河动床模型相似律进行设计，水平比尺、垂直比尺分别为800及80，原型河床比降J=0.000185，选取容重适中的拟焦沙作为模型沙，以便满足水沙运动与河床变形相似条件。试验在平主槽流量2500 m3/s下测得主槽宽度Bc=500 m，主槽平均水深Hc=2.16 m，平均流速Vc=2.31m/s，由式（6）反求出nc=0.01，与原型糙率接近；当流量增加到3300 m3/s时嫩滩上水平均约0.3 m，水面宽B=2500 m，平均水深H=0.73 m，平均流速V=1.80m/s，再由式（6）直接反求出糙率n=0.0061，显然同实际不符，说明对于宽滩窄槽的复式河段，不能按照上述方法计算河流的综合糙率。

20世纪下半叶人们发现水力因子具有不确定性的时候，已注意到糙率的随机性质，鉴于对水流摩阻行为进行精确的确定性描述超出了数学所能实现的范畴，文献[26]将糙率视为随机变量，通过分析水流与河床泥沙运动之间相互影响机理，建立了合理模拟糙率随机变化的概率模型。

1987年秦荣昱[37]以Keulegan的阻力公式为基础，分析了包括黄河等方面实测点据，反求综合糙率Ks，指出在一般情况下Ks＞D65，在有沙波的情况下，Ks比D65大1～3个数量级，还引出了沙波的虚拟当量糙度概念，给出的计算方法颇有意义。1988年Peterson A W与Peterson A E.[38]利用大量实测资料分别回归了弗劳德数Fr（=V（/gH）0.5）小于0.5和Fr＞0.6两个条件下V 同R、J、D的经验关系式，在D对V影响方面两式定性相反，尤其第一式D 竟跟V成正比。1990年王士强以Engelund 研究[36]为基础，通过水槽试验和各种影响因素分析，结合众多实测资料多次回归适线，提出的阻力计算方法包括众多计算式[39]，可适用于不同水力条件下低能态、过渡态及高能态区域，体现了水流阻力问题的复杂性和研究的系统性。

鉴于对数流速分布公式推导中引用了忽略黏性剪切力的假定，使公式在边壁附近不适用，故依此为基础研究水流阻力时就必须进行处理，人为性与运用难度增大，学术性削弱。为在黄河等河流数学模型开发中解决糙率模拟问题[40-41]，不靠直接事先赋值再根据计算结果的需要加以人为调整，赵连军和张红武1995年[42]采用能克服对数流速公式边壁附近计算不合理缺陷同时能反映含沙量影响的张红武流速分布公式[25]，引入摩阻厚度δ*后沿垂线积分，推导出在黄河上应用广泛的如下动床糙率公式：

式中：cn为涡团参数，对于清水取cn= 0.15，对于挟沙水流，可用下式计算：

式中：体积计含沙量Sv=S/γs；γs为泥沙的容重。

对于粗糙区，δ*取当量粗糙度Δ，对于光滑区，δ*取黏性底层厚度δ，并给出过渡区摩阻厚度计算式，从而使上式具有较高的学术价值。但对于天然冲积河流，当量粗糙率Δ 受黄河沙波产生、长消影响较大，受1961年张瑞瑾选用Fr作为研究沙波运动的参变量[43]的启发，并通过动床模型试验按比尺换算成的原型资料，建立摩阻厚度δ*同Fr及D50的关系，即[41-42]：

上式首次在阻力研究中突出了Fr的影响作用，从定量上描述了黄河下游河道当Fr=0.5时沙波消失使δ*最小（床面为动平床），以及Fr＞0.5后因逆行沙波的逐渐出现使δ*复而增加的规律。2001年笔者与彭杨在利用数学模型研究三峡水库与泾河东庄水库泥沙冲淤问题时，为尽量使Fr较小时摩阻厚度δ*计算值较为合理，将当年张红武引入的摩阻厚度δ*的最大值限制为0.1 h，即δ*＞0.1 h后均取δ*=0.1 h。此外，还发现在利用挟沙水流资料建立式（15）时，Fr 已不自觉地反映了水流条件对含沙量的影响，故式（13）计算糙率时可直接取涡团参数cn= 0.15。

1999年Wu和Wang 将Fr、沙粒切应力τ′、临界切应力τc同摩阻参数A 相联系，采用三次多项式拟合水槽及野外实测点据，得到了动床糙率公式[44]：

式中：T=τ′/τc50，τc50为床沙中径对应的临界切应力，为总切应力，n′为沙粒阻力对应的糙率，由前人一系列经验式计算[41]。

2007年邓安军等[45]根据黄河干支流实测资料点绘了n～Fr关系图，并直接以指数形式给出如下拟合关系式：

同时根据王桂仙、Guy及Simons 等学者水槽试验资料分析n～Fr关系，看出“Fr增加到某一临界值后出现n随Fr增大而增大的‘向上翘’现象”，其认识颇有意义。此外，还通过黄河干支流资料检验，认为上式“计算结果与目前在黄河干支流上运用较多的张红武公式比较，两公式计算准确程度基本一样，但该公式形式比张红武公式简便得多”。不过，邓安军采用式（13）计算时，所取“cn=0.95”，比涡团参数cn的正常取值偏大过多，可能在公式其他方面的计算有误。

张罗号[24]2012年系统研究了黄河河槽糙率异常的原因，并给出的修正公式等解决途径颇有意义。2016年Ma和Huang[46]采用黄河下游水文站资料，回归了3个水流阻力经验公式。第1个公式仅考虑水深和比降的影响，第2个公式增加了河宽的影响，第3个公式在第1个公式基础上增加了河宽和含沙量的影响，其研究颇有意义，但验证发现3个公式精度偏低，只能用来定性分析河床特征。2017年麻妍妍和夏军强采用大量资料验证表明，式（13）对天然河道阻力精度最高，并以黄河下游实测资料为依据，建立了糙率与Fr的对数型经验公式[47]：

上式形式简便，具有便于应用的优点。只是文中对于实测数据筛选按照的原则存在问题，如“比降应在（1.5～4）×10-4之间”，岂不知黄河下游孙口、艾山、泺口、利津水文站对应的河床比降与测量的水面比降基本上在（1～1.5）×10-4之间，且花园口、夹河滩、高村这3个处于游荡型河段的水文站，所测量的水面比降很多也在（1～1.5）×10-4之间，导致“剔除有问题的资料”后反而使验证数据代表性受到影响。实际上，大部分学者对阻力公式检验选择的江河资料都缺乏代表性。例如，黄才安等选取国内外河流资料进行验证，选用了我国黄河、长江两条河流的资料[48]。尤其黄河资料，J、D 不少数据突破了正常范围，且Fr=0.015～0.232 而缺乏的Fr=0.25～0.65 常见的资料，长江资料也有类似缺陷（如Fr=0.007～0.026 而缺乏的Fr=0.03～0.3 常见的资料）。

综上所述，现有计算方法存在的问题，一是引入过多系数指数，使理论性与可适用性明显降低。大多在有理论基础的公式处理中，引入较多系数和指数，因事先假定考虑不周而发现资料规律不好后，再视这些系数和指数为某些变量的函数，出现更多系数指数，经验性与运用难度随之变强。例如，上述有几家成果即是在因变量和变量关系式中，将一些无量纲组合参数包括其中，出现不少需水槽资料率定的系数指数，因不符合实际而又假定系数指数同某些变量呈函数形式，利用资料再对新的系数指数加以率定，致使形式繁杂。

二是按水槽试验数据率定系数与指数的公式，因水槽水深很小、弗劳德数较大，同河流实际差距太大。掌握模型相似理论的人知道，水槽试验跟天然原型很难做到严格相似，现实中不易找到相互存在的比尺换算关系，两方面资料直接点绘一起很难拟合，毕竟水槽与水槽、原型与原型的资料自身还难以统一，即使资料放在图上拟合，其曲线也会歪曲实际，故文献[43]认为，不能“轻率地直接运用这些公式去处理实际问题”。此外，进行铺沙的水槽试验难度大、观测精度低，常给出同认识相悖的表现。例如，上述有位学者的水槽试验流量较小时，受尾门调节的比降小至万分之0.1～0.2，Fr很小，测验精度和是否均匀流都难保证，给出了n 跟Fr或V成正比的反常点群关系。加上国外资料包括了大量非冲积河流观测资料，其河床状态对水流阻力的影响规律不正常，也会导致公式误差很大。为此，不少人验证公式时不直接比较n、λ，而是以h（V或C 等）的计算值与实测值加以比较[44，48]。尽管如此，计算h或V时常需试算，某些组次会出现不收敛而被排除[48]，且易出现其他麻烦，如Peterson 等[38]的方法，不同Fr条件利用不同公式试算V 实际很困难，何况他们在Fr=0.5～0.6之间还缺乏计算式。

三是不少阻力计算方法是已知J、V、D 等条件，通过试算出h后再求n或λ，方法不实用，有的公式甚至连h 也是已知的，更显不出必要性。对于天然河流情况，V、h、D 测验精度较高，是易知的条件，最难已知且精度最难保证的是能坡J。毕竟水面比降观测记录中常出现反常现象，黄河、长江水文站曾因此而取消相关测验项目，且目前很大部分河流水文站一直无该测验内容。况且水力能坡J同水面比降、河床比降并不相等[8]，将水面比降作为能坡J属于无奈之举。假如各因子都齐全，直接由曼宁公式推求n 就行了，而不必再绕圈采用Niazkar M[49]之类公式计算n，精度差且概念也有问题。

四是直接通过水沙因子资料回归的经验公式局限性大。以指数型拟合的n～Fr关系式为例，在宽浅断面（H ≈R）可将曼宁公式（6）表示为n 同R1/6J1/2成正比、同Fr 成反比的表达形式，由于验证的n数据本身是按曼宁公式反算的，沙质河流R1/6J1/2变化一般不大，回归公式因受曼宁公式主导，自然会同实测数据符合较好。再将式（17a）表达为：V=0.000731n-1.363(gH)0.5，与形式大不相同的式（6）联解，得（取g=9.8m/s2）：

显然，该式同曼宁公式差距很大，不会符合河流实际。由此看出，现有糙率计算方法还存在明显缺陷，以至于现在工程界水流计算与科研界河流数学模型确定糙率，常常采用经验赋值的办法，影响了学科发展水平。

应该承认，在曼宁公式本身就是纯经验公式的前提下，且糙率测验资料精度实难保证的客观状况下，给出理论性的糙率公式很不现实。按满足工程应用需要，最迫切的是要建立符合精度较高、适用范围广、结构简便三要求的糙率公式。为此，本团队近些年在前人研究基础上，以同原型存在着相似关系的黄河下游河段模型试验观测为关键技术依据，集中对糙率计算方法进行系统攻关，从不同角度都取得了进展。本文先介绍考虑主要影响因素建立的冲积河流糙率计算公式。

3 本文糙率公式的建立与验证

由对数流速沿水深分布公式对流速求导得流速梯度[50]：

式中：z为水深坐标；κ为浑水卡门常数，与式（14）表达的涡团参数cn的关系为：κ=0.375cn，则可表示为[50]：

式中κ0为清水卡门常数，一般取0.4。

于是，在冲积河流一般是挟沙水流的条件下，由上述流速梯度表达式与式（14）看出，κ或cn是反映悬移质含沙量S 引起流速分布变化的特征参数[50]，能体现S 对水流内部能耗的影响。在式（13）至式（15）中，通过cn可反映S 对水流内部能耗的影响，通过D50和Fr可反映水力泥沙因子变化及天然河床上附加形体兴衰对水流摩阻特性的影响。之所以该公式对黄河糙率计算精度最高[47]，其原因是抓住了沙粒肤面摩阻因子D501/6、相对粗糙度D50/h、S、Fr 这4个影响冲积河流糙率的主要因素。其中式（13）与式（14）体现了公式的理论基础，式（15）体现了公式的经验特色。若想提高对其他河流验证精度，以相应资料率定式（15）中系数指数即可实现。

沙波引起的河床糙率是冲积河流糙率的重要组成部分，对河流力学不少课题的研究都有密切关系。但由于沙波产生、发展、消衰表现过于复杂，且难以确切观测，一直难有实质性突破。1925年Exner[51]对单宽输沙率加以简化处理后，对水流、泥沙的连续方程联立求解，得到了沙纹形成、发展的数学表达式。即使以此研究计起，至今也接近一个世纪，但在沙波运动方面的研究仍没取得定量上能供工程应用的进展。在今后相当长时期内，人们还不可能真正搞清冲积河流沙波兴衰过程的机理及沙粒阻力和沙波阻力之间定量关系的客观条件下，当下应着重于建立反映最主要影响因素的动床糙率计算公式。

从实用角度看，最需要满足的就是精度高、适用强、结构形式简单三个条件。为此，根据文献[40]反映的上述4个影响糙率的因素基础上，在沙质河床泥沙起动流速（或临界起动剪切力）变化相对不大条件下，为回避因引入Shields 数、沙粒剪切雷诺数等参数而需知J的麻烦，必须抓住Fr 这个同其他判别沙波运动状态的参数有关联性与可替代性的首要影响参数。1961年Simons 等[52]的室内试验，给出了不同Fr条件下床面出现静平床、沙纹、沙波发展、消衰至动平床等全过程，且上述基于动床比尺模型试验的式（15），也能描述黄河原型沙波产生、长消过程对摩阻厚度δ*的影响规律，说明如此做法是有试验根据的，且同张瑞瑾等[43]以Fr作为研究沙波运动参变量的技术路线一致。

由式（18）可知，对数型关系式难免出现“Fr极小时n会极大、Fr很大时n会极小”甚至偏离原有物理含义的缺陷；采用指数型公式如式（17）除存在类似问题外，还会同曼宁公式产生一定矛盾。封丘河段对于黄河下游具有代表性，原型床沙中径D50≈0.07 mm，将该河段模型测量值按比尺换算成原型资料后发现：n随Fr增大而减小，Fr很小时n会较大；Fr=0.5时糙率n=0.01，河床形态接近动平床，同上述式（15）1995年依据的黄河下游动床模型试验与Simons 等[6，52]的早期试验结果接近（后者动平床的Fr=0.55～67＞0.5，原因是试验采用粒径为0.28 mm的粗沙，泥沙抗拒运动的能力更大，从而沙波发展相对滞后）。为此，参考张罗号等[53]运用河床自动调整作用原理给出细沙河流床沙中径公式的经验，引入a、b两个待定系数，分母以a+bFr的形式来适应动床模型试验反映的n～Fr规律，将公式形式表达为：

式中：a、b为待定系数；n0为参考糙率，可视为动平床状态下的糙率。

为适应上式中bFr 比a 小一个数量级时n会有较大变化，根据国内外糙率表[6-8，16-17]，结合黄河糙率表现[24]，冲积河流糙率最大值一般为糙率较小值（约为动平床糙率）的10倍左右，可取a=0.1（当bFr小至可忽略不计时，n ≈10n0）。由Fr=0.5、n=0.01，可确定待定系数b ≈1.85之后，参考糙率n0≈0.01，则可将糙率简便公式表示为：

利用黄河下游、黄河上游（宁夏、内蒙古河段）、渭河下游、长江干支流（其中汉江下游缺少的D50资料按文献[28]方法近似弥补）及Brownlie 整理的国外河流[54]（包括MIS、ACP、COL、RGR、RIO、CHO、CHP、CHI、RED、MOR、AMC、RGC 等12条河，资料依次为165、140、101、85、38、33、33、32、30、23、10、8组）共2334组国内外河流实测资料进行验证（见表1，Fr=0.018～0.73），表明该公式跟天然河流的实测资料较为符合，尽管公式结构简单，但计算曲线也较好的位于实测点群之中（见图1）。

图1 公式（21）计算曲线与Fr～n点群的关系

在上述基础上，还可合理引入沙粒阻力及含沙量影响因子，提高式（21）的计算精度和适用范围。首先考虑到沙粒糙率是冲积河流动床糙率的重要组成部分，可将式（10）表示的沙粒糙率n同沙粒阻力因子D501/6呈正比关系式（为便于使用，代表粒径取中值粒径），置于式（21）的分子上，为保证公式对于黄河下游计算精度更高一些，按床沙常见的中径D50=0.1 mm计算沙粒糙率n=0.011，需要在分子上引入0.9的系数，以便跟式（21）的分子0.01 相等，则可将考虑沙粒肤面摩阻影响的动床糙率公式表示为：

式中，A50为D50对应的摩阻参数，黄河干支流参照张有龄研究[31]，取A50=19，国外河流与长江参照Strick⁃ler[30]研究并适当考虑粒径较大范围变化的影响，暂取A50=22.5（0.93+300D50），其中D50的单位按m 计。

表1 不同河流实测资料及糙率简便公式验证结果

其次，可通过引入浑水卡门常数来体现含沙量对糙率的影响。鉴于式（21）在上述确定参考糙率的过程中，所依据的都是挟沙水流资料，Fr 已在相当程度上反映了水流条件对含沙量的影响，只考虑n 同（κ / κ0）1/5呈正比即可，故在式（21）分子上引入该影响因素后，进一步能将同时考虑沙粒肤面摩阻与含沙量影响的冲积河流动床糙率计算公式表示为：

利用上述资料对式（23）与式（22）的检验结果也列入表1中，看出考虑沙粒肤面摩阻与含沙量影响后，公式精度进一步得到提高。至于式（23）与式（22）相比，其验证精度提高有限，后者具有结构相对简单的优点，故在糙率测验精度实难保证的状况下，除含沙量较大的黄河下游尽量采用式（23）外，一般情况下均可推荐式（22）作为动床糙率计算公式。

图2 总体资料图示法验证结果

表2a 各糙率公式验证结果（2334组总体数据）

表2b 各糙率公式验证结果（1322组黄河数据）

表2c 各糙率公式验证结果（698组国外河流数据）

图3 黄河资料图示法验证结果

利用这些资料对各家糙率公式进行验证比较结果见表2，一般情况下都有较好的计算精度。图1实际是对糙率简便公式的验证图，表明按图示法对公式的检验具有直观效果。图2与图3分别给出利用总体资料与黄河资料对各家公式的验证结果，能够进一步看出本文公式适用范围相对更大，且点群能更好地分布在45度线两侧。而其他公式还存在这样或那样的问题，如Wu和Wang公式过于复杂，且对于验证数据存在1.03%的剔除率；邓安军等公式与麻妍妍等公式在n＞0.03后计算值偏小，点群主要分布在45度线以下；赵连军及张红武公式用于黄河精度较高，但用于其他河流精度较低。

4 结论

（1）本文在系统回顾前人明渠流速计算公式研究成就的艰辛过程中，以历史文献研究进展为依据，指出较多教科书与现在文献存在的不确切描述这处，以便在水力学及河流学科教学与研究中减少误解，并将河流糙率n定义为“表征河床湿周形状规则性、粗糙程度、河底形态与水沙情况对水流阻滞作用的综合性影响系数”。

（2）在述评糙率计算方法优点与缺点后，指出了现有公式存在的四方面问题。认为在曼宁公式本身就是纯经验公式和糙率测验精度实难保证的状况下，同时在今后相当长时期内人们尚难科学描述沙波兴衰过程的机理及沙粒阻力和沙波阻力之间定量关系的客观条件下，还不可能建立理论性强的糙率公式，最迫切的是建立符合精度较高、适用范围广、结构简便三要求的糙率计算公式。

（3）在沙质河床泥沙起动流速（或临界起动剪切力）变化相对不大条件下，为回避因引入Shields数、沙粒剪切雷诺数等参数而需知能坡J的麻烦，建立公式需抓住Fr 这个同其他判别沙波运动状态的参数有关联性与可替代性的首要影响参数，笔者的黄河模型试验和Simons 等在不同Fr条件下出现沙波兴衰过程的试验以及张瑞瑾以Fr作为沙波运动参变量的成果，可作为如此处理的基础依据。

（4）构建出可适应模型试验反映的“n随Fr增大而减小，Fr很小时n会较大”变化规律的公式形式后，根据冲积河流“糙率最大值一般为糙率较小值（约为动平床糙率）的10倍”的实际情况与“Fr=0.5时糙率接近动平床糙率”的模型试验结果，可确定出待定系数a=0.1、b=1.85，参考糙率n0=0.01，从而建立了冲积河流动床糙率简便公式（21）。

（5）为提高公式的计算精度和适用范围，将沙粒糙率n同沙粒阻力因子D501/6呈正比关系的表达式置于糙率简便公式分子上，从而给出体现沙粒肤面摩擦影响的动床糙率公式（22）。鉴于在上述公式建立的过程中因依据挟沙水流资料而使Fr 已在相当程度上反映了含沙量的影响，故引入卡门常数适当反映含沙量对水流能耗的影响后，进一步给出了同时考虑沙粒肤面摩阻与含沙量影响的动床糙率计算公式（23）。

（6）利用国内外大量冲积河流实测资料进行的系统验证结果表明：本文建立的糙率简便公式（21）与动床糙率公式（23）或式（22）明显比其他公式更符合实际，均可用于工程计算，尤其按图示法进行的检验，更能直观看出本文公式相对于其他公式的优势。除含沙量较大的黄河下游尽量采用式（23）外，一般情况下均可推荐式（22）作为动床糙率计算公式。