張羅號，張紅武，趙晨蘇

（1.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098；2.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084；3.北京科技大學(xué) 土木與資源工程學(xué)院，北京 100083）

1 問題的提出

隨著河流水沙數(shù)學(xué)模型的不斷發(fā)展，推移質(zhì)級配的研究受到廣泛關(guān)注，但現(xiàn)有研究成果無論在計算精度方面還是在天然河流的適用性上尚不完善，主要原因在于推移質(zhì)級配不僅同河床泥沙組成密切相關(guān)，且在一定水流條件下與懸移質(zhì)泥沙交換頻繁，推移質(zhì)泥沙的運動機理十分復(fù)雜［1］。

沖積河流沙質(zhì)河床中，隨著水流條件增強，河床表面可動泥沙的數(shù)量也逐漸增多，并以推移狀態(tài)運動成為推移質(zhì)泥沙，其中細(xì)顆粒還會進(jìn)一步轉(zhuǎn)化為懸移質(zhì)泥沙，此時推移質(zhì)級配常與床沙較為接近［2］。圖1給出了長江漢口水文站推移質(zhì)泥沙及床沙取樣級配曲線的比較結(jié)果，盡管受采樣方法的差異及顆分精度所限，但由該圖仍可看出，沙質(zhì)河床上推移質(zhì)泥沙級配同床沙組成分布狀況頗為相近。

圖1 沙質(zhì)河床推移質(zhì)泥沙與床沙級配實測數(shù)據(jù)的比較

因此，一些學(xué)者通過修正床沙級配得到確定推移質(zhì)級配的半經(jīng)驗方法。具有代表性的研究成果包括李昌華［3］通過將床沙級配中不動泥沙的百分?jǐn)?shù)扣除并標(biāo)準(zhǔn)化后得到的推移質(zhì)級配公式，被稱為“最大粒徑法”；以及張紅武［4］采用推移質(zhì)泥沙平均粒徑與床沙平均粒徑之比修正床沙級配曲線，得到一條與床沙級配曲線平行的推移質(zhì)級配曲線，被稱為“平均粒徑類比法”。這類方法所建立的模式直觀、計算較為簡便，但其缺點是基本沒有考慮水流泥沙運動同推移質(zhì)級配變化的內(nèi)在聯(lián)系。

Gessler［5］根據(jù)統(tǒng)計理論對不同粒徑級下的泥沙的起動概率進(jìn)行了研究，得到推移質(zhì)級配計算式為：

式中：P0i為床沙中第i粒徑組所占沙重百分?jǐn)?shù)；Fi為推移質(zhì)中粒徑小于di的沙重百分?jǐn)?shù)；qdi為粒徑為di泥沙的不動概率；τ0為河床床面水流拖曳力；σ為拖曳力脈動值的標(biāo)準(zhǔn)差；τc為第i粒徑組的泥沙起動臨界拖曳力；x=τ0/τc為積分變量。

上式中泥沙起動條件由均勻沙Shieleds曲線確定［6］，為將該方法延用到非均勻沙中，后來學(xué)者通過研究得到了類似公式，如董永華［7］考慮了非均勻沙起動特點；張啟衛(wèi)［8］假設(shè)了床沙轉(zhuǎn)化為推移質(zhì)泥沙的條件為作用于泥沙顆粒的瞬時底速大于臨界起動流速。此類方法存在的最大局限是必須已知床沙級配曲線，而水文站每年床沙實測資料十分有限，且大部分水文站沒有該項測量內(nèi)容，無法取得床沙與推移質(zhì)泥沙級配同步資料，很難適應(yīng)不斷變化的推移質(zhì)級配的確定。

若直接從物理概念考慮，推移質(zhì)級配與不同分組粒徑的輸沙數(shù)量相關(guān)。其中研究成果包括，樂培久［9］通過建立非均勻沙輸沙率公式得到推移質(zhì)級配計算公式：

式中：gbi為di粒徑分組下的均勻沙推移質(zhì)輸沙率。但由于推移質(zhì)輸沙率測驗精度低而使現(xiàn)有公式多不適用于典型沙質(zhì)河床［2］，故必然導(dǎo)致此類方法建立的公式精度得不到保證。

目前，對挾沙水流中水流紊動作用推移質(zhì)泥沙的規(guī)律尚無定論，但從推移質(zhì)級配影響因素分析而進(jìn)行的理論探討很多。Diplas［10］通過分析泥沙顆粒所受水流作用力大小與作用時間等非恒定因素，并引入沖量的概念來描述推移質(zhì)泥沙顆粒運動規(guī)律。張紅武［11］根據(jù)隨機理論分析了粗沙河床表面紊動特性和推移質(zhì)泥沙顆粒的關(guān)系，建立了粗顆粒推移質(zhì)泥沙級配公式。張緒進(jìn)［12］、陸永軍［13］在大量實測資料基礎(chǔ)上，對推移質(zhì)級配與床沙組成、水沙條件之間的關(guān)系進(jìn)行了探討。此外，何文杜［14］提出了平衡輸沙條件下推移質(zhì)最大粒徑的確定方法。這些基于水力學(xué)及河流動力學(xué)建立的計算方法在理論與使用上都有一定價值，但多數(shù)公式結(jié)構(gòu)過于復(fù)雜，其中部分參數(shù)還需要通過實測資料率定，計算精度也難以保證，且除張紅武公式外，其余均是把推移質(zhì)級配同床沙組成建立關(guān)系，一般僅適用于粗顆粒河床，尤其目前沙質(zhì)河床推移質(zhì)級配曲線更無理論性強且能夠反映水流運動對推移質(zhì)顆粒組成有直接影響的計算公式。為此，本文在分析近壁紊動源區(qū)內(nèi)泥沙運動臨界圖景的基礎(chǔ)上，對沙質(zhì)推移質(zhì)級配分布的表達(dá)形式進(jìn)行了理論探討。

2 推移質(zhì)級配計算方法

2.1 沙質(zhì)河床近壁泥沙運動的臨界圖形我們知道，被概化為二維均勻流的沖積河流，其紊動源區(qū)主要集中在近底由泥沙顆粒組成的床面附近［15］。因受某種擾動或同近壁區(qū)大流速梯度和強剪力相聯(lián)系的壓力差的作用，不斷產(chǎn)生以高頻率、小尺度紊動為主的紊動渦體，這些渦體逐漸離開河底上升擴(kuò)散至全流區(qū)，從而床面附近即成為被L.Prandtl稱之為“渦體作坊”的近壁紊動源區(qū)［16］。天然沙質(zhì)河床的河流，正是在該區(qū)實現(xiàn)河床泥沙與水流相互作用、相互影響的，導(dǎo)致河床組成的不斷調(diào)整或變化。

泥沙學(xué)者將跳躍作為沙質(zhì)推移質(zhì)運動的普遍形式，亦即沙質(zhì)河床近壁紊動源區(qū)的泥沙一旦起動或被水流帶離床面，隨即可能跳起，以跳躍為主要運動形式，且與懸移質(zhì)泥沙存在著相互交換。于是，沙質(zhì)床面顆粒在近壁運動的臨界圖景可概括為：近底流區(qū)床面顆粒被具有瞬時垂向紊動速度的水流帶離床面而起跳，至最大高度后的回落過程中，或在重力作用下以沉速下落，繼續(xù)留在床面；或被路經(jīng)此地的渦體卷走而離開本床面。

由于垂向紊動速度是一個瞬時量，以往研究從動力學(xué)或運動學(xué)觀點出發(fā)建立平衡方程式都是不嚴(yán)格的［17］。鑒于動量是一個瞬時量，且適用于從微觀角度審視紊動流速場與泥沙顆粒的關(guān)系，故在具有水流同泥沙相互作用過程的任何一個瞬間，用該物理量來分析近底顆粒沉浮的臨界條件是合適的。躍動的泥沙在重力作用下以沉速下落的瞬間，如果在垂向遇到相同動量的紊團(tuán)的對撞，即形成相對的平衡，故維持臨界平衡的條件是具有瞬時垂向紊速vb的紊團(tuán)向上的動量等于具有沉降速度ω的泥沙向下的動量。故沿垂向運用動量對撞平衡原理，可將決定近底顆粒沉浮的臨界條件表示為：

式中：m1、m2分別為水流與泥沙相應(yīng)的質(zhì)量。

式（3）之所以對紊動流速引入絕對值符號，系考慮到泥沙沉速為大于零的數(shù)值（實際上，沉速方向向下，只有方向向上的垂向紊速對泥沙的作用才有效）。在連續(xù)介質(zhì)條件下，向上紊動渦團(tuán)與下沉泥沙顆粒的體積應(yīng)該相同，如果以γ、γs分別代表水與泥沙的容重，m1、m2對應(yīng)的密度即分別為γ/g、γs/g，式（3）可表示為：

一般情況下天然沙容重γs是水流容重γ的2.7倍，因此上式表明，沙質(zhì)河床近底顆粒沉浮的臨界條件是水流的瞬時紊動流速等于2.7倍泥沙在水中的沉速，而非兩者相等，同運用動力學(xué)觀點出發(fā)建立的平衡方程式也有差異。

2.2 推移質(zhì)級配公式由于垂向瞬時紊動流速具有高斯分布性質(zhì)［18］，故可給出其概率分布為：

式中：σvb為垂向紊動強度。一般正態(tài)分布密度的系數(shù)分子為1，上式為2是由于瞬時紊速分布只取垂直向上部分；假定向上、向下兩部分的分布相同，也是因為向上部分對泥沙的作用才有效。

由于沙質(zhì)河床泥沙粒徑一般為0.05～2 mm， 范圍涉及過渡區(qū)和滯流區(qū)，而采用現(xiàn)有沉速公式尚不能直接推求出粒徑的顯式，張羅號等［19］利用量綱和諧原理及前人資料，得到如下包括了粒徑范圍為D=0.006～0.9 mm的沉速公式：

式中：ν為運動黏滯系數(shù)，m2/s。

將式（6）代入式（4），并寫成隨機方程：

根據(jù)概率分布函數(shù)計算公式，首先有：

可求出推移質(zhì)粒徑D的分布密度函數(shù)：

數(shù)學(xué)期望：

均方根為：

確定了上述關(guān)系式后，可研究推移質(zhì)組成分布與垂向紊動流速分布的依存關(guān)系。以小于某粒徑Di的泥沙數(shù)目所占泥沙總數(shù)的百分比表示級配曲線（這種方法對于沙質(zhì)河床，所得結(jié)果與重量百分比法基本接近），從概率論的觀點來看，即為：

將式（9）代入上式，可得：

由于被積函數(shù)式（13）中參數(shù)β及垂向紊動強度σvb均為定值，故可由矩形法或梯形法在（0，Di）區(qū)間內(nèi)進(jìn)行數(shù)值積分計算。對于垂向紊動強度，采用如下計算公式［20］：

式中：u*為摩阻流速，m/s；由公式計算；h為水深，m；J為水面比降；z為以河底作為起始點的水深坐標(biāo)，m；Δ為壁面粗糙度，m。

對于粗糙度Δ，本文根據(jù)前蘇聯(lián)學(xué)者Shevelev［21］的試驗成果及張紅武等近些年的試驗資料［22］，得出糙率n與粗糙度Δ的關(guān)系式，即可用實測資料中易得到的糙率值進(jìn)行計算：

通過資料驗證表明，在天然河流常見水深范圍內(nèi)，可利用式（15）求出各河段糙率對應(yīng)的粗糙度。若水深大于5 m且糙率大于0.016，則由卡門紊流粗糙區(qū)沿程損失系數(shù)公式、達(dá)西-魏斯巴赫公式與謝才-曼寧公式［23］聯(lián)立求解，得到的下列關(guān)系式計算：

式中：R為水力半徑，m。

當(dāng)壁面粗糙度Δ數(shù)值很小時不屬于紊流粗糙區(qū)，故不利用式（16）計算。

綜上所述，在具體計算時取河流近底區(qū)z=3Δ，即可由式（14）—（16）求出不同水流條件下的近底紊動強度σvb，隨后由式（9）得到推移質(zhì)級配分布。

3 公式的檢驗

為驗證本文推移質(zhì)級配分布公式在沖積河流沙質(zhì)河床的適用性，采用黃河上游及下游大量實測資料，對本文公式進(jìn)行系統(tǒng)的計算比較（見圖3—圖4）。

由前文圖1表明，沙質(zhì)河床中推移質(zhì)泥沙級配與河床表層床沙顆粒組成資料較為接近，故將缺少推移質(zhì)級配實測資料的黃河上游主要取樣斷面，利用床沙取樣資料代替。表1列舉出黃河寧蒙河段干流所選取的重要河段河槽取樣斷面位置情況［24］，以及河槽表層及河槽深層泥沙顆粒分析結(jié)果中河床泥沙平均粒徑Dcp與中值粒徑D50。由表2也可看到，本文所選取的黃河寧蒙河段主要斷面的床沙組成分布較為均勻，表明推移質(zhì)泥沙與床沙摻混較頻繁。黃河下游均采用推移質(zhì)級配實測資料，由于缺少同步實測水文資料，故選用對應(yīng)水文站/斷面的月平均水力因子（見表2）。

表1 黃河寧蒙河段實測床沙平均粒徑dcp與中值粒徑d50數(shù)據(jù)比較

表2 黃河實測資料水力因子

此外，利用長江沙質(zhì)河段大量實測資料對本文計算方法也進(jìn)行了驗證見圖5，所取驗證資料的水力因子見表3。

本文建立的推移質(zhì)級配計算方法在檢驗時沒有采取經(jīng)驗假定，但從上述檢驗圖中看出，理論計算曲線同實測資料比較接近，計算得到的推移質(zhì)泥沙中值粒徑D50與各斷面所取資料的平均推移質(zhì)中值粒徑基本相等。

4 結(jié)論

圖2 式（13）與黃河上游實測資料的計算比較

圖3 式（13）與黃河中游與下游實測資料的計算比較

表3 長江實測資料水力因子

圖4 式（13）與長江實測資料的計算比較

本文在總結(jié)前人研究成果的基礎(chǔ)上，從分析沙質(zhì)河床近壁紊動源區(qū)的泥沙運動臨界圖景入手，運用垂向動量平衡原理，給出了沙質(zhì)河床推移質(zhì)沉浮的臨界條件，即水流的垂向瞬時紊動流速等于2.7倍泥沙在水中的沉速。在此基礎(chǔ)上列出垂向瞬時紊動流速的隨機方程，并以垂向瞬時紊速具有高斯分布式為條件進(jìn)行求解，并引入垂向紊動強度與近底粗糙度的計算式，建立了理論性強且又不需要已知特征粒徑的推移質(zhì)級配計算方法。

采用黃河、長江主要沙質(zhì)河段大量推移質(zhì)實測資料檢驗結(jié)果表明，本文從理論上建立的方法同天然河流實際頗為符合，可用于計算沖積河流沙質(zhì)河段的推移質(zhì)級配分布曲線，從而可便于提高沙質(zhì)推移質(zhì)輸沙率計算的精度，對于沖積河流推移質(zhì)泥沙的運動也具有重要的理論意義和實用價值，但對于紊動流速場與泥沙顆粒作用機理的理論探討尚需進(jìn)一步研究。

［1］竇國仁.論泥沙起動流速［J］.水利學(xué)報，1960（4）：46-62.

［2］張羅號.沙質(zhì)河床推移質(zhì)輸沙率計算研究［J］.水利學(xué)報，2017，48（4）：467-472.

［3］李昌華.床面上泥沙繞流上舉力系數(shù)的間接確定［J］.泥沙研究，1984（4）：60-63.

［4］張紅武.沖積床面糙率模擬問題的探討［J］.武漢水利電力學(xué)院學(xué)報，1986（3）：92-99.

［5］GESSLER J.The beginning of bedload movement of mixtures investigated as natural armoring in channels［R］.California：W.M.Keck Lab.of Hyd.and Water Res.，Cal.Inst.Tech.，1967.

［6］GESSLER J.Self-stabilizing tendencies of alluvial channels［J］.Journal of the Waterways Harbors&Coastal Engineering Division，1970，96（2）：235-249.

［7］董永華.非均勻推移質(zhì)級配的實驗研究［J］.人民長江，1989（6）：43-49.

［8］張啟衛(wèi).推移質(zhì)級配的計算方法［J］.泥沙研究，1990（4）：41-48.

［9］樂培九.非均勻沙推移質(zhì)輸沙率的研究［J］.水道港口，1991（1）：1-8.

［10］DIPLAS P，DANCEY C L，Celik A O，et al.The role of impulse on the initiation of particle movement under tur-bulent flow conditions［J］.Science，2008，322（5902）：717-720.

［11］張紅武.粗顆粒推移質(zhì)級配的理論計算［C］//水利水電工程靑年學(xué)術(shù)論文集.北京：中國科學(xué)技術(shù)出版社，1992.

［12］張緒進(jìn)，趙世強，陳遠(yuǎn)信.非均勻推移質(zhì)級配研究［J］.泥沙研究，1990（3）：48-55.

［13］陸永軍，張華慶.非均勻沙推移質(zhì)輸沙率及其級配計算［J］.水動力學(xué)研究與進(jìn)展，1991（4）：96-106.

［14］何文社，方鐸，曹叔尤，等.平衡輸沙條件下推移質(zhì)級配［J］.自然科學(xué)進(jìn)展，2002，12（10）：1113-1116.

［15］張瑞瑾，謝鑒衡，等.河流泥沙動力學(xué)［M］.北京：水利電力出版社.1998.

［16］張羅號.基于渦量-動量傳遞理論的天然河流流速與含沙量垂線分布公式［J］.水利學(xué)報，2014，45（4）：566-573.

［17］侯暉昌.河流動力學(xué)基本問題［M］.北京：科學(xué)出版社.1982.

［18］章梓雄，董曾南.黏性流體力學(xué)［M］.北京：清華大學(xué)出版社，2011.

［19］張羅號，張紅武，張錦方，等.泥石流流速計算與模型設(shè)計方法［J］.人民黃河，2015，37（4）：18-24.

［20］張紅武，江恩惠，等.黃河高含沙洪水模型的相似律［M］.鄭州：河南科學(xué)技術(shù)出版社.1994.

［21］SHEVELEV F A.Investigation of the Basic Hydraulic Regularities of Turbulent Motion in Pipes［M］.Kiev：State Publishing House of Literature on Construction and Architecture，1953.（in Russian）

［22］張紅武，李振山，方紅衛(wèi)，等.寧蒙黃河治理對策研究［R］.北京：清華大學(xué)，2016.

［23］張羅號.明渠水流阻力研究現(xiàn)狀分析［J］.水利學(xué)報，2012，43（10）：1154-1162.

［24］安催花，魯俊，吳海亮，等.黃河寧蒙河段分組泥沙起動（河道沖刷）特性研究［R］.鄭州：黃河勘測規(guī)劃設(shè)計有限公司，2016.