黃河河口段河道糙率計(jì)算方法

王萬(wàn)戰(zhàn) 張世安 Jeffrey A Nittrouer

摘要：為了克服固定糙率法低估小流量水位的缺點(diǎn)，基于黃河口利津站糙率一流速實(shí)測(cè)資料，回歸分析得出了糙率（曼寧系數(shù)）計(jì)算公式及河道參數(shù)a的計(jì)算式。通過(guò)一維數(shù)學(xué)模型模擬發(fā)現(xiàn)，在枯水期、非枯水期參數(shù)a值分別采用枯水期參數(shù)a的平均值、全年平均值時(shí)，得出的動(dòng)態(tài)糙率可較好地模擬小流量和大流量水位。參數(shù)a與主槽內(nèi)嫩灘附近地形有較密切的對(duì)應(yīng)關(guān)系，可采用在嫩灘上、下附近取特征水位，進(jìn)而計(jì)算其對(duì)應(yīng)的水力半徑、比降的方法來(lái)估算非枯水期、枯水期的參數(shù)a值。利津一口門床面形態(tài)在較大流量時(shí)多為動(dòng)平整，對(duì)應(yīng)的糙率較小，在糙率變化區(qū)間的下界0.01附近變化;流量較小時(shí)床面形態(tài)多為沙壟，對(duì)應(yīng)的糙率較大，糙率在接近其上界0.035附近變化。在黃河口河道，影響糙率的主要因素是流速和河道水力半徑，床沙粒徑對(duì)糙率影響較小。

關(guān)鍵詞：小流量;水位;水力半徑;流速;一維水動(dòng)力模型;主槽糙率;黃河口

中圖分類號(hào)，TV856;TV882.1 文獻(xiàn)標(biāo)志碼：A

河流阻力問(wèn)題一直是水利部門最為關(guān)心的問(wèn)題之一。水利部門常用曼寧系數(shù)（n）表示糙率。存在的問(wèn)題是國(guó)內(nèi)外對(duì)n進(jìn)行取值或計(jì)算時(shí)考慮的影響因素差別很大，如河道形態(tài)、河床泥沙級(jí)配、植被等[1-3]。因此，影響糙率n的主要因素是什么？有哪些因素影響其在不同的年份、河段發(fā)生變化？本文以黃河河口河段（利津-口門）一維水動(dòng)力模擬為例，對(duì)糙率的變化規(guī)律進(jìn)行分析研究。

1 研究方法

1.1 基本方程

利用一維水動(dòng)力模型模擬利津-口門河段的水流，其基本方程包括連續(xù)方程和圣維南方程：式中：Q為流量，m³/s;q為側(cè)向入流，m³/s;A為過(guò)水?dāng)嗝婷娣e，m²;H為水位，m;R為水力半徑，m;α為動(dòng)量修正系數(shù);C為謝才系數(shù)，C=;g為重力加速度。

1.2 糙率n的計(jì)算公式

利用利津水文站實(shí)測(cè)資料點(diǎn)繪了主槽糙率與水深h的關(guān)系，利用此關(guān)系進(jìn)行了大量模擬，發(fā)現(xiàn)使用糙率n與h的相關(guān)關(guān)系式時(shí)模擬效果較差。點(diǎn)繪糙率n與主槽平均流速v的關(guān)系（見圖1），可知

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

n的確定不僅與流速v有關(guān)，而且還與水力半徑R和比降J有關(guān)，后者可通過(guò)實(shí)測(cè)大斷面和水位資料獲得，此方法通過(guò)參數(shù)a把R、J與v分開，避免了單純依靠經(jīng)驗(yàn)率定參數(shù)的缺陷。

1.3 模型介紹及參數(shù)確定

本文模擬了1990年和2005年利津一口門段的水流。河道平面形態(tài)、橫斷面、水文站、水位站等見圖2（其中1～40為斷面編號(hào)）。受河道演變影響，每個(gè)斷面里程在不同的年份有所不同，但差別不大。1990年、2005年利津一口門段河長(zhǎng)分別為109.6、107.1km。進(jìn)口和出口邊界條件：利津?yàn)檫M(jìn)口邊界，以其流量過(guò)程作為進(jìn)口邊界條件，圖3為1990年、2005年3-11月流量過(guò)程線，1990年、2005年洪峰流量分別為3510、2950m³/s;出口邊界設(shè)在口門處，用口門處海平面作為邊界條件。

初始河道邊界：采用1990年、2005年汛前河道統(tǒng)測(cè)大斷面作為初始邊界，暫不計(jì)算年內(nèi)河床變化。1990年、2005年模擬河段分別有17、40個(gè)斷面，平均斷面間距分別為6.44、2.55km，斷面位置見圖2。2005年典型斷面分別見圖4、圖5（高程為大沽高程）。

水位驗(yàn)證站：利津、一號(hào)壩、西河口有長(zhǎng)期水位觀測(cè)資料，王家莊、十八公里、丁字路口等有短時(shí)段水位觀測(cè)資料。實(shí)測(cè)水位資料用于水面比降、水力半徑計(jì)算和與模擬水位對(duì)比分析。

糙率分別采用固定糙率法和式（1）計(jì)算值。在使用固定糙率法時(shí)，利津一口門段糙率取0.01。使用式（1）時(shí)，依據(jù)逐日沿程水位站實(shí)測(cè)水位求出斷面的水力半徑（R）和比降（J），進(jìn)而求出參數(shù)a。1990年、2005年利津站參數(shù)a的年內(nèi)變化范圍分別為0.012～0.020、0.012～0.022（見圖6）。進(jìn)行數(shù)值模擬時(shí)，參數(shù)a取固定值，分別采用3種情形：取全年平均值、取枯水期平均值、取非枯水期平均值。不考慮冰凌較多的12月—次年2月，2005年調(diào)水調(diào)沙前的3-6月作為枯水期，此后作為非枯水期。2005年利津斷面參數(shù)a在枯水期和非枯水期的平均值分別為0.0148、0.0161，1990年利津斷面參數(shù)a在枯水期和非枯水期的平均值分別為0.0159、0.0169。用同樣的方法求出其他斷面在枯水期和非枯水期的平均值，見圖7。用全年、枯水期和非枯水期參數(shù)a的平均值分別模擬2005年和1990年的水流過(guò)程。

2 不同參數(shù)取法對(duì)模擬水位的驗(yàn)證與分析

2.1 參數(shù)a取年平均值時(shí)水位驗(yàn)證

使用上述邊界條件、參數(shù)a等，按照模型的穩(wěn)定性要求，選擇模擬時(shí)間步長(zhǎng)（30s），分別模擬1990年和2005年的水流情況（見圖8、圖9，水位為大沽高程），可知：

（1）固定糙率法和基于式（1）的動(dòng)態(tài)糙率法在洪峰期間（僅存在于非枯水期）的模擬水位和實(shí)測(cè)水位基本相同。此期間動(dòng)態(tài)糙率法計(jì)算得出的糙率為0.01，與固定糙率法的取值相同。

（2）采用固定糙率法時(shí)，在洪峰前后的小流量時(shí)期，模擬水位明顯低于實(shí)測(cè)水位，差值最大為1.3m。使用動(dòng)態(tài)糙率法時(shí)小流量水位過(guò)程線與實(shí)測(cè)值分兩種情形：①在非枯水期，小流量過(guò)程模擬水位改善明顯，與實(shí)測(cè)水位過(guò)程基本相符，此時(shí)期糙率明顯大于0.01，但是大多數(shù)時(shí)間小于0.035。②在枯水期（見圖9），基于動(dòng)態(tài)糙率法模擬得出的利津和西河口小流量水位與實(shí)測(cè)值基本接近，而其他斷面例如王家莊、一號(hào)壩、丁字路口的水位明顯高于實(shí)測(cè)值，但是此時(shí)對(duì)應(yīng)的糙率值在大部分時(shí)間已經(jīng)達(dá)到最大（0.035），考慮到參數(shù)a的年平均值大于枯水期平均值，說(shuō)明在枯水期參數(shù)a采用年平均值時(shí)計(jì)算的糙率過(guò)大，參數(shù)a取年均值不合適，應(yīng)嘗試采用較小的a值，例如取枯水期的平均值。

2.2 參數(shù)“取枯水期平均值時(shí)水位驗(yàn)證

參數(shù)a取枯水期平均值時(shí)，按照上述方法進(jìn)行模擬，結(jié)果見圖9，可知：

（1）枯水期。利津站水位變化不大。隨著參數(shù)a變小，除了5月底—6月初糙率明顯減小、水位下降較大外，大部分時(shí)間糙率由原來(lái)的0.035轉(zhuǎn)變?yōu)?.035-0.030，利津站水位稍有下降，整體變化不大，即參數(shù)調(diào)整后利津站計(jì)算水位過(guò)程仍與實(shí)測(cè)值基本相符。但是，西河口站糙率明顯減小，水位有所下降，不過(guò)仍比固定糙率（0.01）的計(jì)算誤差小。

王家莊、一號(hào)壩模擬水位過(guò)程改善明顯。隨著參數(shù)a變小，在枯水期除了5月底—6月初糙率明顯減小、水位下降過(guò)大外，大部分時(shí)間兩站糙率明顯減小，水位明顯降低，與實(shí)測(cè)值基本相符。

丁字路口模擬水位有所改善。隨著參數(shù)a變小，大部分時(shí)間糙率減幅較小，糙率稍微低于臨界值（0.035），水位降幅較小，仍高于實(shí)測(cè)值。5月底—6月初糙率明顯減小，水位降低到與實(shí)測(cè)值基本相符的程度。總體來(lái)講，參數(shù)a采用枯水期平均值，比采用年平均值的模擬效果有所改善。

（2）非枯水期。利津、王家莊、一號(hào)壩、西河口、丁字路口站參數(shù)a變小后，洪峰時(shí)段糙率仍處于臨界值（0.01），計(jì)算水位沒(méi)有變化，仍與實(shí)測(cè)值基本相符。但是，小流量過(guò)程的糙率明顯減小，水位低于實(shí)測(cè)值的幅度更大。

總之，參數(shù)a采用枯水期平均值時(shí)，枯水期水位模擬值與實(shí)測(cè)值基本相符，但是非枯水期小流量水位過(guò)程誤差較大。

2.3 參數(shù)a取非枯水期平均值時(shí)水位驗(yàn)證

參數(shù)a的非枯水期平均值明顯高于年平均值，模擬結(jié)果見圖9，可知：

（1）在非枯水期的洪峰時(shí)段計(jì)算的動(dòng)態(tài)糙率仍然為0.01，因此洪峰時(shí)段水位基本沒(méi)有變化，與實(shí)測(cè)水位基本相符。在非枯水期的其他小流量時(shí)段，糙率明顯增大，甚至達(dá)到臨界值0.035，因此大部分時(shí)段水位升高，明顯高于實(shí)測(cè)值。

（2）在枯水期，由于參數(shù)a取年平均值時(shí)糙率在大部分時(shí)段已經(jīng)達(dá)到了臨界值0.035，因此在枯水期取更大的a時(shí)糙率大部分仍是臨界值0.035，計(jì)算水位與參數(shù)a取年平均值時(shí)的計(jì)算水位相比沒(méi)有明顯變化。

綜上可知，在枯水期，參數(shù)a取枯水期平均值能較好地模擬枯水期的水位過(guò)程;在非枯水期，參數(shù)a取年平均值能較好地模擬非枯水期的洪峰和小流量時(shí)段的水位過(guò)程。

3 參數(shù)a的簡(jiǎn)便取法

為了滿足生產(chǎn)中洪水預(yù)報(bào)的需要，筆者進(jìn)一步研究了參數(shù)a的簡(jiǎn)便取法。把圖6參數(shù)a逐日值與年均值相等或相近的日期挑選出來(lái)，畫出對(duì)應(yīng)的水位線于各自的斷面上（見圖4、圖5），發(fā)現(xiàn)這些水位線平交于斷面主槽“嫩灘（或稱一級(jí)灘地）”偏上附近，即圖4、圖5中特征水位1附近，此線以下為枯水河槽。此嫩灘枯水時(shí)出露、大水時(shí)被淹沒(méi)，無(wú)植被生長(zhǎng)。

用同樣的方法，求出對(duì)應(yīng)于枯水期平均值的水位，畫于各自的斷面上（見圖4、圖5），可見枯水期的特征水位2低于嫩灘。所以取對(duì)應(yīng)于嫩灘上、下的特征水位，計(jì)算對(duì)應(yīng)的水力半徑和比降，就能比較容易地估算出非枯水期、枯水期的參數(shù)a值。

4 糙率的影響因素分析

從式（1）可知，曼寧系數(shù)（rz）的主要影響因素是流速和參數(shù)a。

4.1 水力半徑

參數(shù)a由水力半徑和比降組成。研究發(fā)現(xiàn)，比降隨流量變化較小[4-5]，參數(shù)a的變化主要受水力半徑的影響。從圖7可知，2005年路家莊、清2斷面的參數(shù)a值明顯大于各自上下游附近斷面的a值，兩者之間的河段參數(shù)a值明顯偏小，呈現(xiàn)“兩頭大、中間小”的特征。1990年也有類似特征，但是由于1990年測(cè)驗(yàn)斷面較少，因此參數(shù)a最大值的位置與2005年稍有不同。在同樣的流速條件下，a值的大或小意味著糙率可能較大或較小。

4.2 床面形態(tài)

Simons-Richardson床面形態(tài)判別圖的縱軸為水流功率，其等于剪切力乘以平均流速，橫軸為床沙粒徑[6]，利津一口門床沙粒徑呈現(xiàn)時(shí)空變化特征，其中1990年、2005年利津床沙中值粒徑為0.046～0.094mm。分別計(jì)算1990年和2005年枯水期和非枯水期利津以下各斷面逐日水流功率，根據(jù)相應(yīng)的床沙粒徑，結(jié)合床面形態(tài)判別圖，得出各斷面床面形態(tài)隨時(shí)間的變化圖（見圖8、圖9）。由圖8、圖9可知，在非枯水期洪峰和其他較大流量時(shí)段，床面形態(tài)多為動(dòng)平整，此時(shí)段對(duì)應(yīng)的糙率較小，在0.01附近變化;流量較小時(shí)，床面形態(tài)多為沙壟，偶有沙紋出現(xiàn)，對(duì)應(yīng)的糙率較大，在0.035附近變化。

Simons-Richardson床面形態(tài)判別圖反映出在床沙粒徑為0.1～0.2mm時(shí)，床面形態(tài)主要對(duì)水流功率變化敏感，而對(duì)此范圍的床沙粒徑變化不敏感。因此，在黃河河口，河道糙率與流速、河道水力半徑、泥沙粒徑和河床互相作用，互相影響，其中粒徑大小對(duì)糙率的影響較小。

5 結(jié)論

（1）基于黃河口利津站糙率一流速實(shí)測(cè)資料，回歸得出了流量小于平攤流量時(shí)期的糙率計(jì)算公式n=av^-1，河道參數(shù)a=R^2/3J^1/2。通過(guò)一維數(shù)學(xué)模型模擬發(fā)現(xiàn)，在枯水期、非枯水期參數(shù)a值分別采用枯水期參數(shù)a的平均值、全年平均值時(shí)，得出的動(dòng)態(tài)糙率不僅能較好地模擬洪峰時(shí)段的水位，而且能模擬小流量時(shí)的水位過(guò)程，克服了固定糙率法低估小流量水位過(guò)程的缺點(diǎn)。

（2）河道參數(shù)a與主槽內(nèi)嫩灘附近地形有較密切的關(guān)系，可用于在嫩灘上、下附近分別取其特征水位，進(jìn)而計(jì)算其對(duì)應(yīng)的水力半徑、比降來(lái)估算非枯水期、枯水期的參數(shù)a值。

（3）在本文提出的糙率計(jì)算方法中，糙率主要與流速和水力半徑有關(guān)，更重要的是參數(shù)a把河道水力半徑、比降與流速分開，能反映不同年份河道變化對(duì)糙率的影響。粒徑大小對(duì)糙率的影響較小。

（4）利津一口門床面形態(tài)在較大流量時(shí)多為動(dòng)平整，對(duì)應(yīng)的糙率較小，糙率在其變化區(qū)間的下界0.01附近變化。流量較小時(shí)床面形態(tài)多為沙壟，對(duì)應(yīng)的糙率較大，在上界0.035附近變化。

本文提出的動(dòng)態(tài)糙率計(jì)算公式初步考慮了流速、水力半徑等因素，今后還需要進(jìn)行大量的模擬驗(yàn)證，以確定是否還應(yīng)加入其他因素。

參考文獻(xiàn)：

[1]楊振懷.中國(guó)水利百科全書[M].北京：中國(guó)水利水電出版社，2004：142.

[2]李煒.水力計(jì)算手冊(cè)[M].北京：中國(guó)水利水電出版社，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.