黄河河口段河道糙率计算方法

王万战 张世安 Jeffrey A Nittrouer

摘要：为了克服固定糙率法低估小流量水位的缺点，基于黄河口利津站糙率一流速实测资料，回归分析得出了糙率（曼宁系数）计算公式及河道参数a的计算式。通过一维数学模型模拟发现，在枯水期、非枯水期参数a值分别采用枯水期参数a的平均值、全年平均值时，得出的动态糙率可较好地模拟小流量和大流量水位。参数a与主槽内嫩滩附近地形有较密切的对应关系，可采用在嫩滩上、下附近取特征水位，进而计算其对应的水力半径、比降的方法来估算非枯水期、枯水期的参数a值。利津一口门床面形态在较大流量时多为动平整，对应的糙率较小，在糙率变化区间的下界0.01附近变化;流量较小时床面形态多为沙垄，对应的糙率较大，糙率在接近其上界0.035附近变化。在黄河口河道，影响糙率的主要因素是流速和河道水力半径，床沙粒径对糙率影响较小。

关键词：小流量;水位;水力半径;流速;一维水动力模型;主槽糙率;黄河口

中图分类号，TV856;TV882.1 文献标志码：A

河流阻力问题一直是水利部门最为关心的问题之一。水利部门常用曼宁系数（n）表示糙率。存在的问题是国内外对n进行取值或计算时考虑的影响因素差别很大，如河道形态、河床泥沙级配、植被等[1-3]。因此，影响糙率n的主要因素是什么？有哪些因素影响其在不同的年份、河段发生变化？本文以黄河河口河段（利津-口门）一维水动力模拟为例，对糙率的变化规律进行分析研究。

1 研究方法

1.1 基本方程

利用一维水动力模型模拟利津-口门河段的水流，其基本方程包括连续方程和圣维南方程：式中：Q为流量，m³/s;q为侧向入流，m³/s;A为过水断面面积，m²;H为水位，m;R为水力半径，m;α为动量修正系数;C为谢才系数，C=;g为重力加速度。

1.2 糙率n的计算公式

利用利津水文站实测资料点绘了主槽糙率与水深h的关系，利用此关系进行了大量模拟，发现使用糙率n与h的相关关系式时模拟效果较差。点绘糙率n与主槽平均流速v的关系（见图1），可知

n=av^b（0.01≤n<≤0.035）（1）式中：b=-1;a=R²/³J^1/2，J为比降。

n的确定不仅与流速v有关，而且还与水力半径R和比降J有关，后者可通过实测大断面和水位资料获得，此方法通过参数a把R、J与v分开，避免了单纯依靠经验率定参数的缺陷。

1.3 模型介绍及参数确定

本文模拟了1990年和2005年利津一口门段的水流。河道平面形态、横断面、水文站、水位站等见图2（其中1～40为断面编号）。受河道演变影响，每个断面里程在不同的年份有所不同，但差别不大。1990年、2005年利津一口门段河长分别为109.6、107.1km。进口和出口边界条件：利津为进口边界，以其流量过程作为进口边界条件，图3为1990年、2005年3-11月流量过程线，1990年、2005年洪峰流量分别为3510、2950m³/s;出口边界设在口门处，用口门处海平面作为边界条件。

初始河道边界：采用1990年、2005年汛前河道统测大断面作为初始边界，暂不计算年内河床变化。1990年、2005年模拟河段分别有17、40个断面，平均断面间距分别为6.44、2.55km，断面位置见图2。2005年典型断面分别见图4、图5（高程为大沽高程）。

水位验证站：利津、一号坝、西河口有长期水位观测资料，王家庄、十八公里、丁字路口等有短时段水位观测资料。实测水位资料用于水面比降、水力半径计算和与模拟水位对比分析。

糙率分别采用固定糙率法和式（1）计算值。在使用固定糙率法时，利津一口门段糙率取0.01。使用式（1）时，依据逐日沿程水位站实测水位求出断面的水力半径（R）和比降（J），进而求出参数a。1990年、2005年利津站参数a的年内变化范围分别为0.012～0.020、0.012～0.022（见图6）。进行数值模拟时，参数a取固定值，分别采用3种情形：取全年平均值、取枯水期平均值、取非枯水期平均值。不考虑冰凌较多的12月—次年2月，2005年调水调沙前的3-6月作为枯水期，此后作为非枯水期。2005年利津断面参数a在枯水期和非枯水期的平均值分别为0.0148、0.0161，1990年利津断面参数a在枯水期和非枯水期的平均值分别为0.0159、0.0169。用同样的方法求出其他断面在枯水期和非枯水期的平均值，见图7。用全年、枯水期和非枯水期参数a的平均值分别模拟2005年和1990年的水流过程。

2 不同参数取法对模拟水位的验证与分析

2.1 参数a取年平均值时水位验证

使用上述边界条件、参数a等，按照模型的稳定性要求，选择模拟时间步长（30s），分别模拟1990年和2005年的水流情况（见图8、图9，水位为大沽高程），可知：

（1）固定糙率法和基于式（1）的动态糙率法在洪峰期间（仅存在于非枯水期）的模拟水位和实测水位基本相同。此期间动态糙率法计算得出的糙率为0.01，与固定糙率法的取值相同。

（2）采用固定糙率法时，在洪峰前后的小流量时期，模拟水位明显低于实测水位，差值最大为1.3m。使用动态糙率法时小流量水位过程线与实测值分两种情形：①在非枯水期，小流量过程模拟水位改善明显，与实测水位过程基本相符，此时期糙率明显大于0.01，但是大多数时间小于0.035。②在枯水期（见图9），基于动态糙率法模拟得出的利津和西河口小流量水位与实测值基本接近，而其他断面例如王家庄、一号坝、丁字路口的水位明显高于实测值，但是此时对应的糙率值在大部分时间已经达到最大（0.035），考虑到参数a的年平均值大于枯水期平均值，說明在枯水期参数a采用年平均值时计算的糙率过大，参数a取年均值不合适，应尝试采用较小的a值，例如取枯水期的平均值。

2.2 参数“取枯水期平均值时水位验证

参数a取枯水期平均值时，按照上述方法进行模拟，结果见图9，可知：

（1）枯水期。利津站水位变化不大。随着参数a变小，除了5月底—6月初糙率明显减小、水位下降较大外，大部分时间糙率由原来的0.035转变为0.035-0.030，利津站水位稍有下降，整体变化不大，即参数调整后利津站计算水位过程仍与实测值基本相符。但是，西河口站糙率明显减小，水位有所下降，不过仍比固定糙率（0.01）的计算误差小。

王家庄、一号坝模拟水位过程改善明显。随着参数a变小，在枯水期除了5月底—6月初糙率明显减小、水位下降过大外，大部分时间两站糙率明显减小，水位明显降低，与实测值基本相符。

丁字路口模拟水位有所改善。随着参数a变小，大部分时间糙率减幅较小，糙率稍微低于临界值（0.035），水位降幅较小，仍高于实测值。5月底—6月初糙率明显减小，水位降低到与实测值基本相符的程度。总体来讲，参数a采用枯水期平均值，比采用年平均值的模拟效果有所改善。

（2）非枯水期。利津、王家庄、一号坝、西河口、丁字路口站参数a变小后，洪峰时段糙率仍处于临界值（0.01），计算水位没有变化，仍与实测值基本相符。但是，小流量过程的糙率明显减小，水位低于实测值的幅度更大。

总之，参数a采用枯水期平均值时，枯水期水位模拟值与实测值基本相符，但是非枯水期小流量水位过程误差较大。

2.3 参数a取非枯水期平均值时水位验证

参数a的非枯水期平均值明显高于年平均值，模拟结果见图9，可知：

（1）在非枯水期的洪峰时段计算的动态糙率仍然为0.01，因此洪峰时段水位基本没有变化，与实测水位基本相符。在非枯水期的其他小流量时段，糙率明显增大，甚至达到临界值0.035，因此大部分时段水位升高，明显高于实测值。

（2）在枯水期，由于参数a取年平均值时糙率在大部分时段已经达到了临界值0.035，因此在枯水期取更大的a时糙率大部分仍是临界值0.035，计算水位与参数a取年平均值时的计算水位相比没有明显变化。

综上可知，在枯水期，参数a取枯水期平均值能较好地模拟枯水期的水位过程;在非枯水期，参数a取年平均值能较好地模拟非枯水期的洪峰和小流量时段的水位过程。

3 参数a的简便取法

为了满足生产中洪水预报的需要，笔者进一步研究了参数a的简便取法。把图6参数a逐日值与年均值相等或相近的日期挑选出来，画出对应的水位线于各自的断面上（见图4、图5），发现这些水位线平交于断面主槽“嫩滩（或称一级滩地）”偏上附近，即图4、图5中特征水位1附近，此线以下为枯水河槽。此嫩滩枯水时出露、大水时被淹没，无植被生长。

用同样的方法，求出对应于枯水期平均值的水位，画于各自的断面上（见图4、图5），可见枯水期的特征水位2低于嫩滩。所以取对应于嫩滩上、下的特征水位，计算对应的水力半径和比降，就能比较容易地估算出非枯水期、枯水期的参数a值。

4 糙率的影响因素分析

从式（1）可知，曼宁系数（rz）的主要影响因素是流速和参数a。

4.1 水力半径

参数a由水力半径和比降组成。研究发现，比降随流量变化较小[4-5]，参数a的变化主要受水力半径的影响。从图7可知，2005年路家庄、清2断面的参数a值明显大于各自上下游附近断面的a值，两者之间的河段参数a值明显偏小，呈现“两头大、中间小”的特征。1990年也有类似特征，但是由于1990年测验断面较少，因此参数a最大值的位置与2005年稍有不同。在同样的流速条件下，a值的大或小意味着糙率可能较大或较小。

4.2 床面形态

Simons-Richardson床面形态判别图的纵轴为水流功率，其等于剪切力乘以平均流速，横轴为床沙粒径[6]，利津一口门床沙粒径呈现时空变化特征，其中1990年、2005年利津床沙中值粒径为0.046～0.094mm。分别计算1990年和2005年枯水期和非枯水期利津以下各断面逐日水流功率，根据相应的床沙粒径，结合床面形态判别图，得出各断面床面形态随时间的变化图（见图8、图9）。由图8、图9可知，在非枯水期洪峰和其他较大流量时段，床面形态多为动平整，此时段对应的糙率较小，在0.01附近变化;流量较小时，床面形态多为沙垄，偶有沙纹出现，对应的糙率较大，在0.035附近变化。

Simons-Richardson床面形态判别图反映出在床沙粒径为0.1～0.2mm时，床面形态主要对水流功率變化敏感，而对此范围的床沙粒径变化不敏感。因此，在黄河河口，河道糙率与流速、河道水力半径、泥沙粒径和河床互相作用，互相影响，其中粒径大小对糙率的影响较小。

5 结论

（1）基于黄河口利津站糙率一流速实测资料，回归得出了流量小于平摊流量时期的糙率计算公式n=av^-1，河道参数a=R^2/3J^1/2。通过一维数学模型模拟发现，在枯水期、非枯水期参数a值分别采用枯水期参数a的平均值、全年平均值时，得出的动态糙率不仅能较好地模拟洪峰时段的水位，而且能模拟小流量时的水位过程，克服了固定糙率法低估小流量水位过程的缺点。

（2）河道参数a与主槽内嫩滩附近地形有较密切的关系，可用于在嫩滩上、下附近分别取其特征水位，进而计算其对应的水力半径、比降来估算非枯水期、枯水期的参数a值。

（3）在本文提出的糙率计算方法中，糙率主要与流速和水力半径有关，更重要的是参数a把河道水力半径、比降与流速分开，能反映不同年份河道变化对糙率的影响。粒径大小对糙率的影响较小。

（4）利津一口门床面形态在较大流量时多为动平整，对应的糙率较小，糙率在其变化区间的下界0.01附近变化。流量较小时床面形态多为沙垄，对应的糙率较大，在上界0.035附近变化。

本文提出的动态糙率计算公式初步考虑了流速、水力半径等因素，今后还需要进行大量的模拟验证，以确定是否还应加入其他因素。

参考文献：

[1]杨振怀.中国水利百科全书[M].北京：中国水利水电出版社，2004：142.

[2]李炜.水力计算手册[M].北京：中国水利水电出版社，2006：45.

[3]JONATHAN M N.Roughness Characteristics of New ZealandRivers[J].Journal of Hydraulic Engineering，2001，127（6）：526-527.

[4]ROSGEN D.Applied River Morphology[M].2nd ed.Minne-apolis：Printed Media Companies，1996：8-20.

[5]LEOPOLD B L，WOLMAN G M，MILLER P J.FluvialProcesses in Geomorphology（Dover Earth Science）[M].[S.l.]：[s.n.]，1995：522.

[6]RICHARDSON E V，SIMONS D B.Resistance to Flow inAlluvial Channels[J].Journal of the Hydraulics Division，1966，86：73-99.