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

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

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

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

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

1 研究方法

1.1 基本方程

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

1.2 糙率n的計算公式

利用利津水文站實測資料點繪了主槽糙率與水深h的關(guān)系，利用此關(guān)系進(jìn)行了大量模擬，發(fā)現(xiàn)使用糙率n與h的相關(guān)關(guā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)，后者可通過實測大斷面和水位資料獲得，此方法通過參數(shù)a把R、J與v分開，避免了單純依靠經(jīng)驗率定參數(shù)的缺陷。

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

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

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

水位驗證站：利津、一號壩、西河口有長期水位觀測資料，王家莊、十八公里、丁字路口等有短時段水位觀測資料。實測水位資料用于水面比降、水力半徑計算和與模擬水位對比分析。

糙率分別采用固定糙率法和式（1）計算值。在使用固定糙率法時，利津一口門段糙率取0.01。使用式（1）時，依據(jù)逐日沿程水位站實測水位求出斷面的水力半徑（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ù)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年的水流過程。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

用同樣的方法，求出對應(yīng)于枯水期平均值的水位，畫于各自的斷面上（見圖4、圖5），可見枯水期的特征水位2低于嫩灘。所以取對應(yīng)于嫩灘上、下的特征水位，計算對應(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年測驗斷面較少，因此參數(shù)a最大值的位置與2005年稍有不同。在同樣的流速條件下，a值的大或小意味著糙率可能較大或較小。

4.2 床面形態(tài)

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

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

5 結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

[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.