黃河下游河道準(zhǔn)二維泥沙數(shù)學(xué)模型研究

1、    黃河下游河道準(zhǔn)二維泥沙數(shù)學(xué)模型研究            作者：張紅武 趙連軍 王光時(shí)間：2007-11-25 12:07:00                     摘要：本文首先簡(jiǎn)要回顧了以往黃河河道泥沙數(shù)學(xué)模型研究成果，再以

2、水流運(yùn)動(dòng)方程及經(jīng)過(guò)作者修正的泥沙運(yùn)動(dòng)方程為基礎(chǔ)，同時(shí)引入與實(shí)測(cè)資料相符合的水流挾沙力、動(dòng)床阻力、泥沙級(jí)配等計(jì)算公式作為補(bǔ)充方程，構(gòu)造出黃河下游河道準(zhǔn)二維泥沙數(shù)學(xué)模型。然后，采用1986年11月1996年10月這10年長(zhǎng)系列實(shí)測(cè)資料，開(kāi)展了驗(yàn)證計(jì)算。其結(jié)果表明，該模型不僅能計(jì)算黃河下游河道一般洪水引起的河床沖淤變形，還能成功地模擬出大沙年下游處于強(qiáng)烈淤積時(shí)的規(guī)律。 關(guān)鍵詞：黃河下游河道 準(zhǔn)二維泥沙 數(shù)學(xué)模型  1 黃河河道數(shù)學(xué)模型研究的簡(jiǎn)要回顧黃河數(shù)學(xué)模型的研究，是與流域規(guī)劃、工程建設(shè)和管理運(yùn)用等生產(chǎn)緊密結(jié)合的。早在1955年編制黃河綜合利用規(guī)劃技術(shù)經(jīng)濟(jì)報(bào)告中，就曾在黃河三門峽水庫(kù)規(guī)劃

3、階段用初級(jí)的一維恒定平衡輸沙模型對(duì)水庫(kù)淤積和下游河道的河床沖刷變形進(jìn)行計(jì)算1。當(dāng)時(shí)的計(jì)算結(jié)果認(rèn)為，在桃花峪下游沖刷9年后河床刷深27m，沖刷量及沖刷速度顯然比實(shí)際夸大甚多。三門峽水庫(kù)投入運(yùn)用后，庫(kù)區(qū)嚴(yán)重淤積，并迅速向上游延伸，與原來(lái)的計(jì)算結(jié)果截然不同。為研究改建方案，通過(guò)實(shí)測(cè)資料建立了許多泥沙沖淤量與水文要素之間的經(jīng)驗(yàn)關(guān)系式。麥喬威、趙業(yè)安、潘賢娣等學(xué)者在分析下游河道沖淤變化及挾沙能力變化規(guī)律的基礎(chǔ)上，提出了三門峽水庫(kù)下游河床沖淤計(jì)算的方法1。在1975年1977年進(jìn)行治黃規(guī)劃過(guò)程中，為研究各種規(guī)劃方案對(duì)黃河下游的減淤作用，黃河水利委員會(huì)及其有關(guān)科研、設(shè)計(jì)部門的科技人員，使用兩種河道沖淤計(jì)算的

4、經(jīng)驗(yàn)?zāi)Ｐ烷_(kāi)展了大量的計(jì)算工作，麥喬威、李保如兩位學(xué)者進(jìn)行了總結(jié)。其中，在規(guī)劃方案的初步比較階段，使用的是李保如提出的框算模型，只考慮影響輸沙能力的最主要因素（來(lái)水、來(lái)沙）。對(duì)于問(wèn)題比較復(fù)雜和初選的方案，則使用較詳細(xì)的經(jīng)驗(yàn)?zāi)Ｐ?。20世紀(jì)80年代，劉月蘭、韓少發(fā)、吳知等學(xué)者在上述經(jīng)驗(yàn)?zāi)Ｐ偷幕A(chǔ)上，根據(jù)不同河段汛期 、非汛期輸沙經(jīng)驗(yàn)公式，引入水流連續(xù)方程、動(dòng)量方程及泥沙平衡方程式，建立了一套頗具特色的黃河下游河道沖淤計(jì)算模型3。該模型在黃河規(guī)劃中得到應(yīng)用，只是在第三次加固大堤設(shè)計(jì)中，預(yù)計(jì)的下游淤積量與實(shí)際相差甚大。1987年1990年，謝鑒衡、韋直林、魏良琰同劉月蘭、張啟衛(wèi)等，開(kāi)展了黃河下游洪水演

5、進(jìn)預(yù)報(bào)水動(dòng)力學(xué)數(shù)學(xué)模型研究。特別是1992年之后，對(duì)國(guó)家“八五”重點(diǎn)科技攻關(guān)項(xiàng)目研究的全面啟動(dòng)，使黃河數(shù)學(xué)模型研究有了進(jìn)一步發(fā)展的機(jī)遇，在錢意穎、張啟衛(wèi)主持下，開(kāi)展了“黃河泥沙沖淤數(shù)學(xué)模型的應(yīng)用”專題的研究4，在此期間，黃科院、清華大學(xué)、武漢水利電力大學(xué)及中國(guó)水科院的有關(guān)學(xué)者，對(duì)黃河數(shù)學(xué)模型的研究都有較大貢獻(xiàn)。附帶指出，黃委會(huì)有關(guān)專家曾對(duì)美國(guó)最通用的泥沙模型HEC6模型，用黃河資料進(jìn)行了檢驗(yàn)，認(rèn)為這類模型直接用于黃河必須做很多的工作。另外，美國(guó)學(xué)者楊志達(dá)博士等1986年開(kāi)發(fā)了GSTARS模型，參照我國(guó)經(jīng)驗(yàn)進(jìn)行修正后，對(duì)黃河楊集至孫口河段進(jìn)行了驗(yàn)證計(jì)算。以謝鑒衡為首的專家組的驗(yàn)收意見(jiàn)指出：“本模

6、型要應(yīng)用于包括游蕩性河段在內(nèi)的黃河下游，還需進(jìn)一步做更多的研究工作”。此后，世界銀行又委托沈?qū)W文對(duì)該成果進(jìn)行評(píng)價(jià)，提交的書(shū)面意見(jiàn)認(rèn)為：“模型可以用于黃河某些不平衡輸沙問(wèn)題不是十分嚴(yán)重的順直河段，但橫向輸沙很麻煩且需要豐富的經(jīng)驗(yàn)”顯然，海外數(shù)學(xué)模型同黃河實(shí)際尚存在著很大的距離。在黃河平面二維泥沙數(shù)學(xué)模型的研究方面，韋直林5采用有限元法對(duì)黃河下游夾河灘至高村河段的1982年的洪水演進(jìn)過(guò)程進(jìn)行了模擬：程曉陶、劉樹(shù)坤6以二維恒定流與不平衡輸沙理論為基礎(chǔ)，將顯式有限差分法與有限體積法相結(jié)合，并引入張紅武水流挾沙力公式計(jì)算泥沙沖淤變化，研制出夾河灘至高村河段灘區(qū)水沙運(yùn)行數(shù)學(xué)模型。結(jié)合工程實(shí)際，張世奇、楊國(guó)

7、錄等也對(duì)黃河京廣鐵橋以上河段的沖淤規(guī)律進(jìn)行了數(shù)值計(jì)算。與此同時(shí)，張紅武、江恩惠、趙連軍、李煒、鄭邦民、楊明、鐘德鈺、張凌武等學(xué)者聯(lián)合攻關(guān)，開(kāi)展了黃河下游不同河段的平面二維數(shù)學(xué)模型研究。例如，張紅武、李東風(fēng)、許雨新7等應(yīng)用黃河泥沙運(yùn)動(dòng)基本理論研究的最新成果，采用有限元法建立了濟(jì)南河段平面二維水流泥沙數(shù)學(xué)模型。從河床沖淤形態(tài)、流量水位過(guò)程、全域流速場(chǎng)等方面驗(yàn)證結(jié)果可以看出，該數(shù)學(xué)模型計(jì)算結(jié)果與實(shí)測(cè)資料符合較好。應(yīng)該指出，由于黃河問(wèn)題的復(fù)雜性；現(xiàn)有的數(shù)學(xué)模型對(duì)所存在的參數(shù)，多進(jìn)行經(jīng)驗(yàn)處理或用實(shí)測(cè)資料率定。而黃河水沙條件及河床邊界條件如果發(fā)生很大變化，那些經(jīng)驗(yàn)性較強(qiáng)的參數(shù)的處理，往往會(huì)影響新條件下的計(jì)

8、算結(jié)果，因而不能用來(lái)預(yù)測(cè)黃河未來(lái)沖淤趨勢(shì)。另一方面，水流挾沙力及河床糙率是黃河泥沙數(shù)學(xué)模型應(yīng)考慮的兩個(gè)重要問(wèn)題，不少模型中沒(méi)有解決這兩個(gè)問(wèn)題。因此，即使在驗(yàn)證時(shí)差別不大，但在重要的特殊洪水期的計(jì)算結(jié)果，與原型卻差異甚多，大大影響了數(shù)學(xué)模型在黃河防汛規(guī)劃中的使用價(jià)值。再者，恢復(fù)飽和系數(shù)的確定也是黃河泥沙數(shù)學(xué)模型的關(guān)鍵和難點(diǎn)所在，理論值應(yīng)大于1，而一般的黃河模型中的取值都遠(yuǎn)小于1，其處理值得進(jìn)一步商榷。為此謝鑒衡指出：“現(xiàn)階段 對(duì)黃河數(shù)學(xué)模型上的恢復(fù)飽和系數(shù)，無(wú)論理論認(rèn)識(shí)或?qū)嶋H經(jīng)驗(yàn)都甚感不足，有必要在這兩方面做進(jìn)一步的工作?！蔽覀冋J(rèn)為，該問(wèn)題主要是在假定非平衡狀態(tài)下的恢復(fù)飽和系數(shù)與平衡狀態(tài)下相等之

9、后出現(xiàn)的，存在著理論上的缺陷，因此，無(wú)論對(duì)該系數(shù)如何進(jìn)行改進(jìn)，都難以從根本上提高泥沙數(shù)學(xué)模型的理論水平。作者針對(duì)上述存在的問(wèn)題，長(zhǎng)期開(kāi)展了黃河數(shù)學(xué)模型的研制工作，本文介紹黃河下游準(zhǔn)二維泥沙數(shù)學(xué)模型的研究結(jié)果。2 準(zhǔn)二維泥沙數(shù)學(xué)模型的建立2.1 基本方程 本模型計(jì)算除選用描述河道水流運(yùn)動(dòng)的基本方程外，還改正了在泥沙數(shù)學(xué)模型中占據(jù)重要位置的河床變形方程及泥沙連續(xù)方程式8,9，即：河床變形方程泥沙連續(xù)方程式(1)(2)中：i為斷面號(hào)；j為子斷面號(hào)，河床高程最低的子斷面號(hào)j為1，最高的取j為m；A為過(guò)水面積；t為時(shí)間；x為沿流程坐標(biāo)；Z為水位；為泥沙渾水沉速；S為水流含沙量；S為水流挾沙力；0為河床淤

10、積物干容重；bij為子斷面寬度；Zbij為子斷面平均河床高程；為平衡含沙量分布系數(shù)，由下式計(jì)算10： 其中：式中：為卡門常數(shù)；cn為渦團(tuán)參數(shù)(cn=0.375)；C為謝才常數(shù)；u為摩阻流速；且有： 運(yùn)用量綱分析及相似轉(zhuǎn)化原理，并根據(jù)動(dòng)床河工模型試驗(yàn)結(jié)果，可求得Ks的表達(dá)式為： 對(duì)于非飽和系數(shù)fs，通過(guò)歸納分析，再借助試驗(yàn)資料，得出：2.2 水流挾沙力及糙率計(jì)算 為保證數(shù)學(xué)模型適用于黃河下游高含沙洪水，提高多沙河流輸沙計(jì)算的精度，水流挾沙力計(jì)算采用黃河上常用的張紅武水流挾沙力公式10 式中： 式中：sk為第k組粒徑泥沙在渾水中群體沉速；ok為第k組粒徑沉速；D50為床沙中值粒徑；d50為懸沙中值

11、粒徑(mm)；psk為第k組粒徑泥沙的重量百分?jǐn)?shù)；N為非均勻沙的分組數(shù)。經(jīng)舒安平、江恩惠、朱太順及陳雪峰、陳立、李義天11等學(xué)者，通過(guò)黃河、長(zhǎng)江等河流大量實(shí)測(cè)資料系統(tǒng)檢驗(yàn)結(jié)果表明，式(8)是現(xiàn)有公式中最適用于黃河等天然河流的水流挾沙力公式，目前已得到廣泛應(yīng)用。而且教科書(shū)中也被推薦為適用于高含沙水流的公式12。本模型采用作者的動(dòng)床糙率公式開(kāi)展黃河洪水水位模擬計(jì)算，可較好地描述水力泥沙因子的變化對(duì)摩阻特性的影響，并能反映天然河道中各種附加糙率的影響。該糙率公式為：13 式中：為摩阻厚度，采用下游資料，得出與床沙中    徑D50及佛汝德數(shù)Fr的關(guān)系式為： 2

12、.3 懸移質(zhì)泥沙與床沙級(jí)配的計(jì)算 從懸浮于水體中泥沙顆粒的受力分析入手，經(jīng)理論推導(dǎo)和經(jīng)驗(yàn)處理可建立懸移質(zhì)泥沙級(jí)配及粒徑二階圓心距的計(jì)算公式14 上述3式中，P(di)為小于某粒徑的泥沙所占的百分?jǐn)?shù)；d50為懸沙中值粒徑；dcp為懸沙平均粒徑；m為指數(shù)；d為懸沙粒徑的二階圓心距，其數(shù)值可表征出泥沙組成的非均勻程度。由式(15)、(16)可知，對(duì)于變量d50、dcp、d，已知任意兩個(gè)，即可確定第三個(gè)變量的值，懸沙組成即可由式(14)求解出來(lái)。同時(shí)，通過(guò)黃河下游 水文站實(shí)測(cè)資料回歸分析，建立了懸沙中值粒徑d50和平均粒徑dcp之間的如下關(guān)系式：式中：vm為渾水運(yùn)動(dòng)粘滯系數(shù)。對(duì)于天然細(xì)沙河流，底沙隨時(shí)

13、可能跳起，同懸沙的交換幾率很大。我們以躍動(dòng)模式分析床沙起動(dòng)時(shí)的受力入手，推導(dǎo)出與懸沙相類似的級(jí)配及粒徑二階圓心距公式： 式中：P(Di)為床沙中小于某粒徑的泥沙所占的百分?jǐn)?shù)；Dcp為床沙平均粒徑；m為指數(shù)；D為床沙粒徑的二階圓心距。從懸沙與床沙交換過(guò)程中任一粒徑組的沙量平衡入手，經(jīng)理論推導(dǎo)建立了如下懸沙與床沙平均粒徑與二階圓心距變化規(guī)律的表達(dá)式：式(21)式(24)中：dL、L分別為側(cè)向入流泥沙的平均粒徑及粒徑的二階圓心距，P為床沙混合層厚度，dc、c為沖淤物平均粒徑與粒徑的二階圓心距，dc與P的具體計(jì)算方法參考文獻(xiàn)14。通過(guò)對(duì)黃河下游實(shí)測(cè)資料回歸得出如下懸沙平均粒徑沿橫向分布公式： 式中：d

14、cpi為斷面平均懸沙平均粒徑，dcpij為斷面上任一點(diǎn)懸沙平均粒徑。C2為斷面形態(tài)系數(shù)，由沙量守衡可求得 2.4 懸移質(zhì)含沙量橫向分布計(jì)算 通過(guò)對(duì)黃河大量實(shí)測(cè)資料分析發(fā)現(xiàn)，含沙量的橫向分布規(guī)律不僅與水力因子、含沙量大小有關(guān)，還與懸沙組成密切相關(guān)。懸沙粒徑越細(xì)，含沙量的橫向分布越均勻，為此除引入含沙量因子外，還引入懸浮指標(biāo)/ku*，建立了含沙量橫向分布公式，即 式中：Vi、Vij分別為斷面平均及任意一點(diǎn)的流速；Hi、Hij分別為斷面平均及任意點(diǎn)的水深；u為斷面平均摩阻流速；C1為1左右的斷面形態(tài)系數(shù)，由沙量守恒可求得：式中：qij為斷面任一點(diǎn)單寬流量，y為橫向坐標(biāo)；a、b為斷面河寬兩端點(diǎn)起點(diǎn)距。

15、2.5 河槽在沖淤過(guò)程中河寬變化規(guī)律的模擬 本數(shù)學(xué)模型所采用的模擬河槽寬度變化過(guò)程的方法，是根據(jù)計(jì)算河段應(yīng)受特定的河相關(guān)系均衡調(diào)整原理進(jìn)行的。模擬河段內(nèi)各橫斷面寬度修正值，除受到給定河岸條件的限制外(例如受到抗沖河岸、山嘴及河道整治工程等的限制)，還要隨著河流造床過(guò)程而自動(dòng)調(diào)整。本模型選用張紅武河床綜合穩(wěn)定性指標(biāo)9作為河相關(guān)系均衡調(diào)整準(zhǔn)則計(jì)算，即：式中：B為河寬；i為河床縱比降。作為初步探討，將河床的均衡關(guān)系式運(yùn)用于河床斷面的自動(dòng)調(diào)整過(guò)程，對(duì)于河寬(包括子斷面寬度)的調(diào)整以式(29)的約束關(guān)系確定其變化值。3 驗(yàn)證計(jì)算結(jié)果 該數(shù)模建立之后，主要針對(duì)黃河下游諸如“82.8”、“77.7”、“92

16、.8”、“96.8”等不同類型洪水開(kāi)展了驗(yàn)證計(jì)算。為進(jìn)一步檢驗(yàn)該數(shù)學(xué)模型對(duì)不同水沙系列年的適用性，特選用三門峽水庫(kù)蓄水運(yùn)用及滯洪運(yùn)用初期1960年9月1965年6月的水沙系列和近期的1986年11月1996年10月水沙系列以及相應(yīng)的沖淤變形資料，對(duì)模型進(jìn)行檢驗(yàn)。前者見(jiàn)文獻(xiàn)15，本文只介紹后者的驗(yàn)證計(jì)算結(jié)果。1986年1996年，黃河下游總來(lái)水量為2979.6億m3，總沙量為85.742億t，其中汛期來(lái)水量為1404.8億m3，占總來(lái)水量的47.1%，來(lái)沙量為81.7億t，占總來(lái)沙量的95.3%。特別是1988年、1992年、1994年、1996年，這4年為豐沙年，汛期小浪底水文站平均含沙量均超

17、過(guò)80kg/m3，汛期來(lái)沙量為48.6億t，占總來(lái)沙量的56.6%。由于水沙搭配不合理，造成了下游河道嚴(yán)重淤積，10年內(nèi)下游河道共淤積25.0億t，其中鐵謝至高村河段淤積19.2億t,占總淤積量的76.8%，高村至艾山河段淤積3.6億t，占總淤積量的14.4%，艾山至利津河段淤積2.2億t，占總淤積量的8.8%。從淤積的年際年內(nèi)分配上看，河道主要在汛期淤積，共淤積31.4億t，而非汛期河道基本處于沖刷狀態(tài)，共沖刷6.4億t。其中1988年、1992年、1994年、1996年四年汛期淤積最為嚴(yán)重，共淤積22.4億t，占總淤積量的89.5%。從泥沙的淤積部位看，大量泥沙主要淤積在主槽內(nèi)，引起河槽高

18、程明顯抬升，平灘流量大幅度下降。如“92.8”、“96.8”兩場(chǎng)洪水，花園口站洪峰流量分別為6020m3/s、7600m3/s，接近于該站洪峰的多年平均值，僅相當(dāng)于中常洪水，花園口站卻先后出現(xiàn)有水文記載以來(lái)的最高水位，洪水漫灘范圍及傳播時(shí)間也明顯增大，洪水造成的經(jīng)濟(jì)損失巨大。采用小浪底站實(shí)測(cè)水沙過(guò)程及懸沙組成作為進(jìn)口邊界條件。伊洛河、沁河入黃水沙及泥沙級(jí)配過(guò)程分別采用黑石關(guān)、小董(或武陟)站實(shí)測(cè)資料，沿黃引水引沙量采用實(shí)測(cè)資料沿程分配給出，作為已知的分流條件。初始地形采用黃河下游1986年汛后92個(gè)實(shí)測(cè)大斷面控制，并將每一個(gè)大斷面概化為1215個(gè)子斷面。初始床沙組成根據(jù)花園口、夾河灘、高村、孫

19、口、艾山、濼口、利津等7站的實(shí)測(cè)床沙組成按斷面間距內(nèi)插或外延。出口邊界水位采用相應(yīng)年份 利津站水位流量關(guān)系控制。表1給出了鐵謝至利津不同河段不同時(shí)段計(jì)算和實(shí)測(cè)的沖淤量。計(jì)算表明，1986年11月1996年10月，鐵謝至高村河段淤積13.588億m3，高村至艾山淤積2.445億m3，艾山至利津淤積1.654億m3，與相應(yīng)的實(shí)測(cè)值13.707億m3、2.581億m3、1.571億m3接近。從上述計(jì)算結(jié)果可以看出，無(wú)論河道處于沖刷、淤積或沖淤交替狀態(tài)，即使是遇到大沙年(包括其中的高含沙洪水)處于強(qiáng)烈淤積時(shí)，該模型均能很好地模擬黃河河道的輸沙特性及河床變形規(guī)律。同時(shí)，該模型對(duì)河道在不同水沙組合條件下發(fā)生的上沖下淤、上淤下沖及沿程沖刷或淤 積等規(guī)律的模擬精度也令人滿意。4 結(jié) 語(yǔ)本文首先回顧了前人關(guān)于黃河河道泥沙數(shù)學(xué)模型研究成果，再以黃河泥沙基礎(chǔ)理論最新成果 為基礎(chǔ)，構(gòu)造出黃河下游河道準(zhǔn)二維泥沙數(shù)學(xué)模型。然后采用原型實(shí)測(cè)的水沙過(guò)程及相應(yīng)的河床變形資料，開(kāi)展了驗(yàn)證計(jì)算。結(jié)果表明：該模型不僅能模擬下游河道長(zhǎng)系列年一般洪水引起的河床沖淤變形，還能成功復(fù)演大沙年下游處于強(qiáng)烈淤積時(shí)的實(shí)際情況。該數(shù)學(xué)模型已在黃河下游防洪及規(guī)劃中應(yīng)用，而且還在小浪底水庫(kù)運(yùn)用