二維水沙數(shù)學(xué)模型的建立和應(yīng)用

1、研 究 生：關(guān)建偉指導(dǎo)教師：李大鳴 教授天津大學(xué)建筑工程學(xué)院二零一零年六月摘摘 要要 二十世紀(jì)70年代以來(lái)，隨著電子計(jì)算機(jī)的普及與計(jì)算機(jī)性能的不斷提高，人們?cè)谘芯亢恿髂嗌硢?wèn)題時(shí)，越來(lái)越多的采用了數(shù)學(xué)模型。它以數(shù)值方法和計(jì)算機(jī)技術(shù)為手段，通過(guò)對(duì)河流的數(shù)值模擬計(jì)算，解決河流工程所關(guān)心的問(wèn)題。對(duì)二維水流泥沙數(shù)學(xué)模型的探索和研究已在水利、環(huán)境、港口航道、防洪等工程領(lǐng)域中取得了一定成果，但水流泥沙運(yùn)動(dòng)機(jī)理仍是工程問(wèn)題中的重要研究課題。 本文結(jié)合河北省水利水電勘察設(shè)計(jì)研究院河北省廟宮水庫(kù)汛期調(diào)度運(yùn)用方式項(xiàng)目，研究了平面二維水流泥沙數(shù)學(xué)模型的理論、方法和實(shí)際應(yīng)用問(wèn)題，主要工作如下： 1.系統(tǒng)地推導(dǎo)了二維非恒

2、定流、二維泥沙連續(xù)方程，應(yīng)用有限體積法思想，采用分類(lèi)簡(jiǎn)化離散方程的方法，降低了計(jì)算難度、節(jié)省了計(jì)算時(shí)間,為建立平面二維水沙數(shù)學(xué)模型奠定了基礎(chǔ)。 2.在已有的泥沙研究成果的基礎(chǔ)上，提出并推導(dǎo)了新的泥沙連續(xù)方程形式，結(jié)合適合廟宮水庫(kù)情況的愛(ài)因斯坦水流挾沙公式，建立了平面二維水沙數(shù)學(xué)模型。 3.根據(jù)離散方程，編制Fortran計(jì)算程序，對(duì)廟宮水庫(kù)庫(kù)區(qū)內(nèi)的水流和泥沙沖淤進(jìn)行計(jì)算并與已有的實(shí)測(cè)資料對(duì)比，吻合效果較好，說(shuō)明模型可以用于該水庫(kù)的河道變形模擬。對(duì)上游來(lái)水、來(lái)沙情況和水庫(kù)排沙控制水位序列13年水沙組合過(guò)程進(jìn)行計(jì)算模擬，給出了全河段的流場(chǎng)分布、含沙量分布、河床沖淤演變分布的計(jì)算模擬成果。 4.通過(guò)

3、綜合對(duì)比，定量分析了不同計(jì)算組合下床面沖淤演變情況，計(jì)算出排沙比,提出了水庫(kù)現(xiàn)狀泄流條件下最優(yōu)清淤方案。 研究成果不僅為理論研究提供了新思路，而且也為實(shí)際工程提供了可行的數(shù)值模擬方法,具有一定的理論創(chuàng)新和應(yīng)用價(jià)值。 關(guān)鍵詞關(guān)鍵詞：平面二維水流泥沙數(shù)學(xué)模型；廟宮水庫(kù)；泥沙連續(xù)方程；平面二維水流泥沙數(shù)學(xué)模型；廟宮水庫(kù)；泥沙連續(xù)方程；沖淤演變沖淤演變目 錄第一章 緒論 1.1泥沙數(shù)學(xué)模型發(fā)展簡(jiǎn)況 1.2水沙數(shù)學(xué)模型基本方程及計(jì)算方法 1.3數(shù)值計(jì)算方法 1.4水沙基本理論存在的問(wèn)題 1.5淤積上延現(xiàn)象第二章 河道復(fù)雜邊界的處理技術(shù) 2.1 斜對(duì)角笛卡爾方法的基本原理 2.2斜對(duì)角笛卡爾方法的優(yōu)點(diǎn) 2

4、.3邊界條件及邊界網(wǎng)格單元類(lèi)型處理 2.4 本章小節(jié)第三章 平面二維水流泥沙數(shù)學(xué)模型的建立 3.1 二維水沙數(shù)學(xué)模型的基本方程 3.2 初值條件及邊界條件的處理 3.3 網(wǎng)格劃分 3.4 水流挾沙能力第四章 平面二維水沙數(shù)學(xué)模型的數(shù)值模擬 4.1 水流控制方程的離散 4.2 泥沙連續(xù)方程的離散 4.3 河床變形方程的離散 4.4 模型求解過(guò)程第五章 水流泥沙數(shù)學(xué)模型驗(yàn)證及應(yīng) 用 5.1 模型的基本特征 5.2 網(wǎng)格劃分 5.3 模型驗(yàn)證 5.4 模型預(yù)報(bào) 5.5 排沙量及沖淤分析第六章 工作與展望 6.1 工作 6.2 展望目前研究河流泥沙及河床演變問(wèn)題，主要有三種手段，即原型資料分析法、物理

5、模型方法、數(shù)學(xué)模型模擬方法。這三種方法各有千秋，視具體工程問(wèn)題，可獨(dú)立運(yùn)用或聯(lián)合運(yùn)用。物理模型方法是將原型的尺寸以及各種水力要素在遵循相似原則的基礎(chǔ)上按一定比尺縮小為模型。物理模型方法具有直觀(guān)、準(zhǔn)確、易于操作等特點(diǎn)，是早期研究河床演變的主要方法。近年來(lái)隨著測(cè)量技術(shù)的進(jìn)步，物理模型的精度得到了很大的提高，但是與數(shù)學(xué)模型相比，物理模型耗資大，對(duì)外界條件要求高，且完成實(shí)驗(yàn)方案的周期長(zhǎng)。河流泥沙數(shù)學(xué)模型能夠模擬水流泥沙運(yùn)動(dòng)過(guò)程，河床沖淤量、河床形態(tài)及河床泥沙組成的變化。它的研究和應(yīng)用可以追溯到上世紀(jì)50年代，由于受到當(dāng)時(shí)計(jì)算工具的限制，這種計(jì)算僅限于一些簡(jiǎn)單的狀況，而且精度較差、工作量大，因此，未能得

6、到普遍推廣。二十世紀(jì)70年代以來(lái)，隨著電子計(jì)算機(jī)的普及與計(jì)算機(jī)性能的不斷提高，這些困難已逐步克服，數(shù)學(xué)模型的優(yōu)越性越來(lái)越突出。經(jīng)過(guò)近二、三十年的發(fā)展，水流泥沙數(shù)學(xué)模型經(jīng)歷了由一維到二維、到二 第一章第一章 緒論緒論1.1 泥沙數(shù)學(xué)模型發(fā)展概括泥沙數(shù)學(xué)模型發(fā)展概括 維嵌套、到三維、到三維嵌套，由原先在概化水文、泥沙及河床條件下數(shù)學(xué)模型的建立、率定和驗(yàn)證，發(fā)展到與物理模型相互配合共同回答工程水流和泥沙問(wèn)題(復(fù)合模型)的階段。1.2 水沙數(shù)學(xué)模型 1.水流控制方程 2. 泥沙控制方程0yxqqZtxy222227 322()xxyyxxxxxxtqqqqqqqqqqZghgnthxhyxhxy222

7、227322()yxyyyyyyyxtq qqqqq qqqqZghgnth xh yyhxy*2222()()()()()()kxkykskskskkkksxsyhSq Sq SSStxyhShSxy *1()Nbs ks kkkkZSSt1.3 數(shù)值計(jì)算方法概述1.3.1 數(shù)值離散方法描述河流水沙運(yùn)動(dòng)的控制方程一般是非線(xiàn)性的，對(duì)于這些控制方程，一般很難求得解析解，只能通過(guò)離散控制方程、求解代數(shù)方程組來(lái)得到近似數(shù)值離散解。數(shù)值計(jì)算方法的精度及速度依賴(lài)于控制方程的離散方法、代數(shù)方程組的求解方法、網(wǎng)格的劃分及邊界條件的處理等。目前常用的數(shù)值計(jì)算方法有：特征線(xiàn)法、有限差分法、有限元法 、有限體積法

8、、 有限分析法等。 1.特征線(xiàn)法特征線(xiàn)法是河流模擬數(shù)值計(jì)算在計(jì)算機(jī)尚未普遍之時(shí)所采用的一種方法，主要采用手工計(jì)算，其最初思路是在平面上繪制特征線(xiàn)。其優(yōu)點(diǎn)是能反映問(wèn)題中信息沿特征傳播的性質(zhì)，算法符合水流運(yùn)動(dòng)的物理機(jī)制，穩(wěn)定性好、計(jì)算精度高，但對(duì)于求解周期短、變化急劇的問(wèn)題計(jì)算效果較差，因此目前很少直接用于數(shù)值計(jì)算。但是其原理仍很重要，經(jīng)常用于作為了解其他數(shù)值方法的基礎(chǔ)。2. 有限差分方法有限差分方法(FDM)是計(jì)算機(jī)數(shù)值模擬最早采用的方法，至今仍被廣泛運(yùn)用，F(xiàn)DM以Taylor級(jí)數(shù)展開(kāi)等方法，把控制方程中的微商用差商代替進(jìn)行離散，從而建立代數(shù)方程組來(lái)求解該方法數(shù)學(xué)概念直觀(guān)、表達(dá)簡(jiǎn)單，其解的存在性

9、、收斂性和穩(wěn)定性早已有較完善的研究成果，是比較成熟的數(shù)值方法，目前應(yīng)用廣泛。 3.有限元方法有限元法是根據(jù)極值原理(變分或加權(quán)余量法)，將問(wèn)題的控制微分方程化為控制所有單元的有限元方程，把總體的極值作為各單元極值之和，即將局部單元總體合成，形成嵌入了指定邊界條件的代數(shù)方程組，求解該方程組就得到各節(jié)點(diǎn)上待求的函數(shù)值。有限元的求解方法常見(jiàn)的有直接法、變分法、加權(quán)余量法和能量平衡法等。其優(yōu)點(diǎn)是網(wǎng)格劃分靈活、擬合復(fù)雜河岸邊界容易、網(wǎng)格節(jié)點(diǎn)可局部加密、穩(wěn)定性好、精度高,適合于幾何、物理?xiàng)l件復(fù)雜的問(wèn)題。對(duì)于隱式FEM，其精度較高，但數(shù)學(xué)推導(dǎo)繁雜、計(jì)算量和儲(chǔ)存量較大，而且在誤差估計(jì)、收斂性和穩(wěn)定性等方面的理

10、論研究與有限差分法相比還顯得不夠成熟和完善。尤其在多維計(jì)算中，由于有限元法貯存量大、會(huì)直接影響計(jì)算速度。 4.有限體積法有限體積法是計(jì)算出通過(guò)每個(gè)控制體邊界沿法向輸入(出)的流量和動(dòng)量通量后，對(duì)每個(gè)控制體分別進(jìn)行水量和動(dòng)量平衡計(jì)算，便得到計(jì)算時(shí)段末各控制體平均水深和流速。因此，有限體積法正是對(duì)于推導(dǎo)原始微分方程所用控制體的回歸。因?yàn)榭刂企w間界面輸運(yùn)的通量，對(duì)相鄰控制體來(lái)說(shuō)大小相等、方向相反，故對(duì)整個(gè)計(jì)算域而言，沿所有內(nèi)部邊界的通量相互抵消。對(duì)由一個(gè)或多個(gè)控制體組成的任意區(qū)域，以至整個(gè)計(jì)算區(qū)域，都嚴(yán)格滿(mǎn)足物理守恒律，不存在守恒誤差，并且能正確計(jì)算間斷。由于采用守恒型的微分方程并對(duì)每一計(jì)算單元進(jìn)行

11、質(zhì)量和動(dòng)量守恒形式的離散，使得微分方程包含的守恒性質(zhì)在每一個(gè)控制容積上都得到滿(mǎn)足，若保持各單元界面兩側(cè)相鄰控制體的計(jì)算輸運(yùn)通量相等，那么整個(gè)計(jì)算區(qū)域上都能保持守恒。在數(shù)值模擬的過(guò)程中，雖然采用的離散求解方法不同，但都有相同的特點(diǎn)，即首先把計(jì)算區(qū)域劃分成許多控制體或網(wǎng)格，然后在這些小控制體上把微分方程離散成代數(shù)方程，再把小控制體上的代數(shù)方程匯合成總體代數(shù)方程組，最后在一定的初始條件下求解此方程組，從而求得計(jì)算區(qū)域內(nèi)各節(jié)點(diǎn)的物理量。所以數(shù)值模擬的正確性和精度取決于網(wǎng)格的劃分、方程的離散、初始值條件、代數(shù)方程組的求解以及所建模型的物理理論依據(jù)是否正確合理等幾個(gè)因素。各種方法均有其自身的優(yōu)點(diǎn)和適應(yīng)性，

12、在實(shí)際計(jì)算時(shí)選擇什么數(shù)值方法應(yīng)根據(jù)所研究問(wèn)題的特點(diǎn)和計(jì)算精度要求，以及研究者的習(xí)慣而定。1.3.2邊界處理技術(shù)1.不規(guī)則邊界的處理在水利工程中，水流的邊界往往是很復(fù)雜的。而在數(shù)學(xué)模型計(jì)算中如何處理不規(guī)則的復(fù)雜邊界將直接影響到整個(gè)模型計(jì)算結(jié)果的可靠性、精度和計(jì)算時(shí)間，是數(shù)值模擬的一大難點(diǎn)，選擇復(fù)雜邊界的處理方法與算法的選擇有著同等重要性。常用的不規(guī)則復(fù)雜邊界處理方法主要有：美國(guó)學(xué)者Thompson16提出的貼體坐標(biāo)法、非結(jié)構(gòu)性網(wǎng)格以及多重網(wǎng)格法。貼體坐標(biāo)法因其能使計(jì)算域邊界上的坐標(biāo)線(xiàn)與邊界線(xiàn)密切貼合，而獲得了廣泛的研究與應(yīng)用。其主要優(yōu)點(diǎn)是用曲線(xiàn)坐標(biāo)逼近復(fù)雜邊界，在邊界模擬上具有較高的精度。但該方

13、法在較復(fù)雜的邊界上不易形成正交貼體網(wǎng)格，如果是非正交曲線(xiàn)坐標(biāo)系，則會(huì)產(chǎn)生大量的交叉導(dǎo)數(shù)項(xiàng)，給求解帶來(lái)很大的困難。非結(jié)構(gòu)性網(wǎng)格使用三角形和四邊形計(jì)算單元，使得在復(fù)雜區(qū)域形成計(jì)算網(wǎng)格更加靈活和方便。但非結(jié)構(gòu)性網(wǎng)格往往要建立一套標(biāo)明網(wǎng)格節(jié)點(diǎn)之間關(guān)系的復(fù)雜系統(tǒng)，因此在編程計(jì)算處理時(shí)比較麻煩，特別是針對(duì)自由面流動(dòng)和干濕邊界變化時(shí)，這種困難更加突出。多重網(wǎng)格法則由于網(wǎng)格系統(tǒng)之間的相對(duì)獨(dú)立性，難以做到對(duì)系統(tǒng)守恒性的保證。另外對(duì)于復(fù)雜邊界的處理問(wèn)題，俞國(guó)良于1987年提出了鏡像法，周建軍于1988年提出了邊界對(duì)稱(chēng)點(diǎn)法。這兩種方法都是根據(jù)對(duì)稱(chēng)原理，假設(shè)流場(chǎng)外有一定的虛擬流動(dòng)，從而導(dǎo)出計(jì)算域外節(jié)點(diǎn)各物理量的計(jì)算式

14、。該方法簡(jiǎn)單易行，提高了不規(guī)則邊界上運(yùn)用差分方法的計(jì)算精度，但由于域內(nèi)外對(duì)稱(chēng)點(diǎn)一般不在節(jié)點(diǎn)上，仍需一定插值過(guò)程，會(huì)帶來(lái)一定誤差。 2.動(dòng)邊界的處理由于河床沖淤變化和水位的波動(dòng)，在河岸及洲灘等區(qū)域都會(huì)遇到動(dòng)邊界問(wèn)題。動(dòng)邊界數(shù)值模擬的困難主要在于沿動(dòng)邊界法向流動(dòng)不同于明渠均勻流，其非恒定和非均勻性強(qiáng)，常用的曼寧摩阻公式等在形式上難以套用，再者是因?yàn)樗詈苄?，?duì)離散格式的要求很高，要求數(shù)值解不產(chǎn)生假振，保證水深總是大于零。常見(jiàn)的處理動(dòng)邊界方法有兩種：一種是對(duì)實(shí)際邊界的數(shù)值模擬，即時(shí)刻追蹤動(dòng)邊界的準(zhǔn)確位置。另一種方法是：把整個(gè)計(jì)算區(qū)域都按照常進(jìn)行計(jì)算，在計(jì)算過(guò)程中取一定的方法對(duì)無(wú)水區(qū)域進(jìn)行處理。但是由

15、于第一種方法在編制程序中處理較復(fù)雜，因此目前經(jīng)常采用的方法是第二種方法，在計(jì)算過(guò)程中不追蹤動(dòng)邊界的準(zhǔn)確位置，把整個(gè)計(jì)算域都參與運(yùn)算，對(duì)于干河床部分采用一些處理技巧進(jìn)行處理，比如何少苓提出的“窄縫法”，其思想是在床面上人為設(shè)置窄縫，由于縫足夠深，總是有水流動(dòng)，而縫很窄，當(dāng)河床露出水面時(shí)由窄縫引起的水量守恒誤差很小。此外，還有周建軍提出的滲透邊界法，程文輝、王船海等運(yùn)用的“凍結(jié)”技巧。1.4水沙基本理論中存在的問(wèn)題1.阻力問(wèn)題 阻力計(jì)算是阻力系數(shù)(糙率及渦流粘性系數(shù))的確定問(wèn)題。處理阻力的方法有兩種，一種是采用水力半徑分割法或能坡分割法劃分阻力單元，分別計(jì)算各單元的阻力，然后計(jì)算總阻力；另一種方法

16、是根據(jù)實(shí)測(cè)資料直接計(jì)算總阻力。以上兩種方法都只能用于一維水流糙率的計(jì)算。多數(shù)的二維數(shù)學(xué)模型均是直接采用一維阻力系數(shù)進(jìn)行計(jì)算，這種處理方式顯然比較粗略。李義天在考察一維情況下斷面綜合糙率沿河寬變化分布的情況下，提出了二維糙率計(jì)算公式。實(shí)際應(yīng)用表明，該方法比直接采用一維方法確定的糙率更合理，能比較好地反映實(shí)際情況，但由于因素比較復(fù)雜，由此得到的推論在定量上尚未能被實(shí)測(cè)資料所證實(shí)。2.水流挾沙力問(wèn)題水流挾沙力的計(jì)算一般有兩種方法，一種是將含沙量沿垂線(xiàn)分布公式積分求垂線(xiàn)平均含沙量，但河底泥沙含沙量很難確定，不同的研究學(xué)者給出了不同的計(jì)算方法，所得的計(jì)算結(jié)果與實(shí)際情況差別較大。另一種是半經(jīng)驗(yàn)半理論方法，

17、這樣的公式很多，公式中包含兩類(lèi)參數(shù)，一類(lèi)是反映水力條件的參數(shù)，如流速水深等；另一類(lèi)是反映泥沙特性的參數(shù)，如沉速、泥沙粒徑等。這種方法是建立在實(shí)測(cè)資料的基礎(chǔ)上得，因而結(jié)果比較接近實(shí)際情況，得到了廣泛的應(yīng)用?，F(xiàn)在挾沙力的研究成果多基于一維問(wèn)題，對(duì)二維問(wèn)題研究的較少，而大多數(shù)的二維問(wèn)題也是直接采用一維挾沙力公式進(jìn)行計(jì)算的，但是李義天通過(guò)經(jīng)驗(yàn)比較，發(fā)現(xiàn)二維挾沙力與一維挾沙力有很大的差異。 3.床沙懸移質(zhì)沿程變化的問(wèn)題目前在數(shù)學(xué)模型計(jì)算中存在的尚未很好地解決的問(wèn)題是床沙質(zhì)及床沙級(jí)配的沿程變化問(wèn)題。對(duì)于沖淤甚小的處于自然狀態(tài)下的河流，在許多情況下床沙及其級(jí)配在沖淤過(guò)程中可以假定基本不變，而沖泄質(zhì)級(jí)配則視來(lái)

18、水來(lái)沙條件而定；對(duì)于沖淤幅度甚大的河流，對(duì)床沙質(zhì)及床沙級(jí)配的沿程變化不能不進(jìn)行考慮?，F(xiàn)在大多數(shù)模型將推移質(zhì)和懸移質(zhì)分開(kāi)，將懸移質(zhì)中床沙質(zhì)與沖泄質(zhì)分開(kāi)的做法基本相同，至于將床沙質(zhì)進(jìn)一步分級(jí)來(lái)求不同粒徑的沖淤變化，一般都比較粗糙。韓其為23提出了一套考慮泥沙級(jí)配變化的計(jì)算模式，其基本思想是不區(qū)分床沙質(zhì)與沖泄質(zhì)，假定水流挾沙力的級(jí)配和實(shí)際輸移的級(jí)配一致，而懸移質(zhì)和床沙質(zhì)級(jí)配在每一時(shí)段內(nèi)的變化都看成本時(shí)段內(nèi)河段沖淤的直接后果，計(jì)算斷面每個(gè)粒徑組的含沙量須反復(fù)試算才能確定，該方法對(duì)泥沙輸移的物理考慮比較細(xì)致但計(jì)算過(guò)于復(fù)雜。1.5淤積上延現(xiàn)象當(dāng)水流進(jìn)入到水庫(kù)回水區(qū)后，由于水流條件的改變，流速減小、泥沙落淤

19、，而落淤的結(jié)果反過(guò)來(lái)又促成水流條件的進(jìn)一步改變，使得回水曲線(xiàn)在淤積區(qū)及其上下游一定范圍內(nèi)普遍抬高，為泥沙繼續(xù)在原淤積區(qū)落淤并同時(shí)向上下游發(fā)展提供條件。這種淤積自回水末端(通常指正常蓄水位回水末端)向上游發(fā)展的現(xiàn)象，稱(chēng)為水庫(kù)淤積的上延現(xiàn)象，也就是通常所說(shuō)的“翹尾巴”現(xiàn)象。按淤積中剖面形態(tài)的類(lèi)型，可以分為三角洲淤積、帶狀淤積、和錐體淤積見(jiàn)。05101520253035010203012345相對(duì)高程(米)距 壩 里 程 (公 里 )1-尾部段；2-頂坡段；3-前坡段；4-異重流淤積段；5-壩前淤積段0204060801001200102030405060相對(duì)高程(米)距壩里程(公里)1-原平均河底

20、高程；2-淤積后平均河底高程12圖 1-1 官?gòu)d水庫(kù)淤積縱剖面圖 1-2 豐滿(mǎn)水庫(kù)淤積縱剖面2130740750760高程(米)距壩里程(公里)211-原河床；2-淤積后灘面圖1-3 陜西黑松林水庫(kù)淤積縱剖面 前的3.74變緩為1.8，淤積末端距大壩13.7km，上延系數(shù)1.55。隨著淤積末端不斷上延，上游河床逐年抬高，距大壩9.6km的四合永公路橋處河床淤高3m多，使河床由原來(lái)的地面以下到現(xiàn)在與兩岸持平，甚至高出兩岸耕地及村基，產(chǎn)生嚴(yán)重浸沒(méi)。浸沒(méi)區(qū)內(nèi)數(shù)千畝耕地鹽漬化、沼澤化，不能正常耕種，近千戶(hù)居民住房地基沉陷、墻垣斷裂甚至倒塌，給群眾生產(chǎn)帶來(lái)嚴(yán)重?fù)p失，以至出現(xiàn)影響社會(huì)穩(wěn)定的問(wèn)題，并嚴(yán)重威脅

21、四合永鎮(zhèn)的安全，受到省、部領(lǐng)導(dǎo)和有關(guān)部門(mén)的高度重視。 為延緩廟宮水庫(kù)淤積浸沒(méi)的發(fā)展，自1987年至1992年連續(xù)6年進(jìn)行了空庫(kù)排沙原型實(shí)驗(yàn)。試驗(yàn)中對(duì)不同水沙條件，空庫(kù)排沙時(shí)間，不同排沙方式等進(jìn)行了全面的試驗(yàn)，總結(jié)分析了水庫(kù)現(xiàn)有泄流排沙設(shè)施條件下的減緩水庫(kù)淤積浸沒(méi)的可能性和運(yùn)用效果。1993年由水利部天津院完成的廟宮水庫(kù)空庫(kù)排沙試驗(yàn)物理模型研究和同年由天津大學(xué)泥沙研究室完成的廟宮水庫(kù)排沙洞前期工程泥沙數(shù)學(xué)模型初步研究，為當(dāng)時(shí)增設(shè)新排沙洞位置、洞徑和底高程的確定提供了理論依據(jù)。 近年來(lái)隨著水環(huán)境關(guān)注意識(shí)的提高和水資源高效利用思路的確立，如何在水庫(kù)保水和水庫(kù)排沙間選擇優(yōu)化的水庫(kù)運(yùn)行模式是需要進(jìn)行研究

22、和解決的重要問(wèn)題。根據(jù)廟宮水庫(kù)設(shè)計(jì)正常蓄水位778m(大沽高程，以下不加說(shuō)明處均采用大沽高程)和庫(kù)區(qū)死水位768m的要求，優(yōu)選即能泄水排沙又能保水興利的調(diào)水調(diào)沙方案是本次研究的主要目的。研究過(guò)程擬采用庫(kù)區(qū)限蓄水位768m、770m、772m、774m、776m、778m六級(jí)水位標(biāo)準(zhǔn)，來(lái)水來(lái)沙過(guò)程采用19581959、19711975、19851989和1990的四個(gè)系列三種組合來(lái)模擬不同年份水沙過(guò)程，在18個(gè)計(jì)算系列中進(jìn)行方案優(yōu)選。第二章 河道復(fù)雜邊界的處理技術(shù)自然界中，水流的邊界往往是十分復(fù)雜的，而在水流數(shù)學(xué)模型計(jì)算中如何處理不規(guī)則的復(fù)雜邊界將直接影響到整個(gè)模型計(jì)算結(jié)果的可靠性、精度和計(jì)算時(shí)

23、間，是數(shù)值模擬的一大難點(diǎn)。復(fù)雜邊界的處理方法與算法的選擇有著同等的重要性，因此在研究江河湖泊的水流數(shù)值計(jì)算過(guò)程中，必須要對(duì)邊界作出一個(gè)適當(dāng)?shù)奶幚矸椒ā?.1 斜對(duì)角笛卡爾方法的基本原理斜對(duì)角笛卡爾方法基本原則就是按就近逼近原則，把距實(shí)際水邊線(xiàn)較近的網(wǎng)格節(jié)點(diǎn)連線(xiàn)以擬合實(shí)際水邊線(xiàn)。生成網(wǎng)格時(shí)，根據(jù)邊界是否固定，可采取不同的判定方法。當(dāng)邊界為固定邊界時(shí)，可把固定邊界離散為若干離散節(jié)點(diǎn)，確定各離散節(jié)點(diǎn)坐標(biāo)，按距離最小原則選定固定邊界上各離散節(jié)點(diǎn)的逼近網(wǎng)格節(jié)點(diǎn)，然后依序連接各網(wǎng)格節(jié)點(diǎn)，即形成擬合邊界；對(duì)于動(dòng)邊界，由于邊界線(xiàn)變化不定，則可借助網(wǎng)格節(jié)點(diǎn)間水深的干濕變化和水深值，參照相應(yīng)網(wǎng)格節(jié)點(diǎn)間的地形高差，

24、即可判定距離實(shí)際水邊線(xiàn)較近的網(wǎng)格節(jié)點(diǎn)，合并縱橫向的擬合網(wǎng)格節(jié)點(diǎn)，即可形成擬合瞬時(shí)水邊線(xiàn)。圖2-1 斜對(duì)角笛卡爾法邊界擬合示例2.2斜對(duì)角笛卡爾方法的優(yōu)點(diǎn)斜對(duì)角笛卡爾方法的優(yōu)點(diǎn) 1.邊界擬合效果良好，更接近實(shí)際邊界。由圖2-2可以看出，計(jì)算水邊線(xiàn)與實(shí)際水邊線(xiàn)擬合良好，其邊界長(zhǎng)度誤差一般小于6%，僅為鋸齒法長(zhǎng)度誤差的，且隨網(wǎng)格間距的縮小而更為精細(xì)。圖2-3為按鋸齒法(也以距離最近為原則)確定的計(jì)算水邊線(xiàn)，對(duì)比圖2-2、圖2-3可知斜對(duì)角笛卡爾法所擬合的水邊線(xiàn)比鋸齒法擬合的水邊線(xiàn)更為貼近實(shí)際水邊線(xiàn)。事實(shí)上，斜對(duì)角法就是在鋸齒法基礎(chǔ)上通過(guò)斜對(duì)角連線(xiàn)添加一個(gè)三角形網(wǎng)格單元，使邊界形態(tài)更接近真實(shí)形態(tài)，且避

25、免了鋸齒狀邊界的產(chǎn)生，使流場(chǎng)計(jì)算更精確。2.處理簡(jiǎn)單，附加計(jì)算量小在矩形網(wǎng)格水域邊界上運(yùn)用斜對(duì)角笛卡爾法，計(jì)算時(shí)，只須對(duì)三角形邊界網(wǎng)格單元進(jìn)行處理，而三角形邊界網(wǎng)格單元若把斜邊中點(diǎn)視為一角點(diǎn)，則仍可視為四邊形網(wǎng)格，可用統(tǒng)一方程求解，與內(nèi)部計(jì)算網(wǎng)格無(wú)異，并不像貼體坐標(biāo)擬合一樣須對(duì)基本方程進(jìn)行坐標(biāo)變換而增添許多交叉導(dǎo)數(shù)項(xiàng)，因此，用斜對(duì)角笛卡爾法進(jìn)行邊界處理，不僅保留了離散簡(jiǎn)便、計(jì)算效率高的優(yōu)點(diǎn)，同時(shí)又使邊界擬合更加精細(xì)，避免了鋸齒狀邊界對(duì)流場(chǎng)計(jì)算的影響。實(shí) 際 邊 界 線(xiàn)斜 對(duì) 角 法 擬 合 線(xiàn)實(shí) 際 邊 界 線(xiàn)鋸 齒 法 擬 合 線(xiàn)圖2-2斜對(duì)角笛卡爾法邊界擬合 圖2-3鋸齒法邊界擬合2.3邊

26、界條件及邊界網(wǎng)格單元類(lèi)型處理根據(jù)斜對(duì)角笛卡爾法的逼近原則，在平面二維計(jì)算中擬合水邊線(xiàn)上形成的計(jì)算邊界網(wǎng)格單元存在八種基本類(lèi)型，見(jiàn)下圖。 圖2-4八種網(wǎng)格單元基本類(lèi)型圖2-5 網(wǎng)格類(lèi)型及差分點(diǎn)分布斜對(duì)角笛卡爾方法不僅可以應(yīng)用于有限差分法中，也可應(yīng)用于有限體積法等其它離散格式中，它對(duì)復(fù)雜邊界條件和復(fù)雜水流條件的適應(yīng)性和計(jì)算精確性將使復(fù)雜邊界的處理變得簡(jiǎn)單高效而且準(zhǔn)確，為復(fù)雜邊界問(wèn)題和近岸工程問(wèn)題的計(jì)算提供了一種新思路。第三章第三章 平面二維水流泥沙數(shù)學(xué)模型的建立平面二維水流泥沙數(shù)學(xué)模型的建立天然河道水沙數(shù)學(xué)運(yùn)動(dòng)一般都屬于三維運(yùn)動(dòng)，運(yùn)動(dòng)要素既沿程變化，又沿水深和河寬發(fā)生變化28。由于三維水流運(yùn)動(dòng)比較

27、復(fù)雜，難于計(jì)算，常將運(yùn)動(dòng)要素沿水深方向平均，把三維問(wèn)題轉(zhuǎn)化為平面二維問(wèn)題。平面二維非恒定流方程，也就是淺水方程，它描述的對(duì)象是符合以下條件的均勻流體：(1) 具有自由表面。(2) 以重力為主要的作用力，以水流與固定邊界之間及水流內(nèi)部的摩阻力為主要耗散力，有時(shí)還存在水面氣壓場(chǎng)風(fēng)應(yīng)力以及地球柯氏力等的作用。(3) 水平流速沿垂線(xiàn)近似均勻分布，不考慮實(shí)際存在的對(duì)數(shù)或指數(shù)等形式的垂線(xiàn)流速分布(4) 水平運(yùn)動(dòng)尺度遠(yuǎn)大于垂直運(yùn)動(dòng)尺度，垂向流速及垂向加速度可以忽略，從而水壓強(qiáng)接近于靜水壓強(qiáng)分布。在實(shí)際問(wèn)題中，通?？梢宰鳛榉呛愣魈幚淼木哂凶杂杀砻娴膶?shí)際水流，有如下情況：(1)水深相對(duì)較淺。可以根據(jù)水深h和波

28、長(zhǎng)L之比區(qū)別，當(dāng)0.4時(shí)稱(chēng)為淺水長(zhǎng)波。(2) 水底坡度較緩。設(shè)底坡傾角為，判斷緩坡的條件是。此時(shí)，可以忽略底坡引起的垂直速度和垂直環(huán)流，也不必考慮垂直加速度及由此產(chǎn)生的動(dòng)水壓力。(3) 水面漸變且坡度較緩。既使?jié)M足以上兩個(gè)條件，在特定的情況下，由于受到天然或人工控制，有時(shí)水面比降仍很大，水深沿程變化很大的急變流，水壓力不能用靜水壓強(qiáng)分布來(lái)近似。(4) 無(wú)明顯垂直環(huán)流。由于外力改變或地質(zhì)環(huán)境的改變?cè)诖怪逼矫鎯?nèi)產(chǎn)生二次流，會(huì)給淺水模型帶來(lái)誤差，需要加以改正。此時(shí)方程不能沿垂向積分，要按不可壓流來(lái)處理。密度沿垂向變化較大，改變了靜壓分布，需要作為密度流處理。因此，通常的淺水流動(dòng)指在重力作用下密度均勻

29、，具有自由表面，流動(dòng)近似水平的長(zhǎng)波傳播現(xiàn)象。 相比一維水沙數(shù)學(xué)模型，平面二維水流泥沙數(shù)學(xué)模型對(duì)沿水深平均的水沙各運(yùn)動(dòng)要素(垂線(xiàn)平均流速、懸移質(zhì)含沙濃度、推移質(zhì)輸沙率)以及河床變形等隨流程沿河寬兩個(gè)方向的變化情況均可以進(jìn)行比較細(xì)致的模擬。1.水流控制方程3.1二維水沙數(shù)學(xué)模型的基本方程0yxqqZtxy222227 322()xxyyxxxxxxtqqqqqqqqqqZghgnthxhyxhxy222227 322()yxyyyyyyyxtqqqqqqqqqqZghgnthxhyyhxy2.泥沙控制方程*2222()()()()()()kxkykskskskkkksxsyhSq Sq SSStx

30、yhShSxy *1()NbsskskkkkZSSt 3.關(guān)于懸沙連續(xù)方程的一些探討泥沙在水流中分布就是懸沙的含沙量，而分布的變化即為懸沙含沙量的變化，它是可以用泥沙連續(xù)方程來(lái)描述的。懸沙含沙量的變化不外乎是由對(duì)流項(xiàng)(上游輸沙)，和本地的紊動(dòng)擴(kuò)散項(xiàng)引起的，由此可以把對(duì)流項(xiàng)，紊動(dòng)項(xiàng)引起的懸沙含沙量變化用兩個(gè)連續(xù)方程來(lái)表示。主要思想如下：由三維泥沙擴(kuò)散方程：其中：是單元體內(nèi)的含沙量，可以這樣認(rèn)為單元體內(nèi)含沙量是由上游來(lái)沙和本地床面紊動(dòng)揚(yáng)沙引起的。我們把上游來(lái)沙影響用來(lái)表示，稱(chēng)為輸沙含沙量。本地床面紊動(dòng)揚(yáng)沙用來(lái)表示，稱(chēng)為揚(yáng)沙含沙量。這樣，可以把單元體總的含沙量用如下式子表示： ()()()xyzSS

31、SSSSSuvtxyzxxyyzz SSS 分別根據(jù)時(shí)均，脈動(dòng)的概念，可以得到兩個(gè)方程： 可見(jiàn)，輸沙含沙量和揚(yáng)沙含沙量都是滿(mǎn)足三維泥沙擴(kuò)散方程的。三維方程在物理意義描述上是十分準(zhǔn)確的，但是在實(shí)際應(yīng)用中，因?yàn)橛?jì)算條件在河道發(fā)生沖淤變化時(shí)，懸移質(zhì)底層及床面揚(yáng)沙含沙量的變化是引起床面沖淤的主要因素，而輸沙含沙量只是上游來(lái)沙在本地的一種擴(kuò)散效應(yīng)，在當(dāng)時(shí)刻沖淤變化過(guò)程中并不會(huì)起到舉足輕重的作用。我們把含沙量區(qū)分為輸沙含沙量與揚(yáng)沙含沙量也是為了在計(jì)算沖淤變化過(guò)程中體現(xiàn)這種關(guān)系，以使結(jié)果更符合實(shí)際情況。()()()xyzSSSSSSSuvtxyzxxyyzz()()()xyzSSSSSSSuvtxyzxxy

32、yzz對(duì)輸沙含沙量方程和揚(yáng)沙含沙量二維化()()xySSShShShuhvhtxyxxyy *()()xySS SSSShShShhuhvtxyxxyy S水 面Z0Z0+haSaSbSAB輸 沙 含 沙 量揚(yáng) 沙 含 沙 量如上所述，懸移質(zhì)中造床質(zhì)所占比例是要受條件限制的，同樣的，床面泥沙能否起動(dòng)，以及起動(dòng)后是以何種形態(tài)運(yùn)動(dòng)，這在模型計(jì)算中，也是要予以考慮的。 Shield(1936)對(duì)各種泥沙顆粒進(jìn)行了臨界起動(dòng)實(shí)驗(yàn)，實(shí)測(cè)得到了無(wú)量綱臨界起動(dòng)剪切應(yīng)力與顆粒雷諾數(shù)的關(guān)系曲線(xiàn)， 泥沙顆粒的懸浮指標(biāo)是通過(guò)懸浮表達(dá)式來(lái)計(jì)算， *Zu (3-111) 式中：k是卡門(mén)系數(shù)，是泥沙顆粒的沉速，*u是摩阻流

33、速。 根據(jù)Einstein的觀(guān)點(diǎn)，對(duì)于給定的水流和泥沙條件，可取5z 作為泥沙是否進(jìn)入懸浮狀態(tài)的臨界判斷值。 222222222222220SSSSSSuvuvuvtxyx ty tx y (3-117) 設(shè)式子(3-108)的解為1sin()kiiiiSAatb xc y，代入式子(3-108)，每一項(xiàng)均可以寫(xiě)為(為簡(jiǎn)單起見(jiàn)去掉下標(biāo)i)： 22222sin()(222)0Aatbxcy ab uc vabuacvbcuv (3-118) 如果式(3-108)解存在，則上式(3-109)要滿(mǎn)足如下條件： 222222220ab uc vabuacvbcuv (3-119) 即 2()0abuc

34、v (3-120) 則 abucv (3-121) 其中解sin()SAatbxcy中，, ,A b c為待定系數(shù)，由初始條件和邊界條件確定，為定性給出解S的變化，給定特定參數(shù)vucbaA,，可以畫(huà)出床面定性變化趨勢(shì)。 式(3-108)是橢圓方程，通解肯定是正弦余弦函數(shù)的線(xiàn)性組合，為了定性顯示解得形式，并取0t 時(shí)刻，取式(3-113)、(3-114)為S的解。 2sin(23 )2cos(2 )Sxyxy (3-122) 2sin(75 )4cos(46 )Sxyxy (3-123) 圖(3-2)、圖(3-3)分別顯示式(3-122)、式(3-123)值的變化以及等值線(xiàn)的分布。為清楚起見(jiàn)，圖

35、(3-2)中Z放大比例為20，, x y空間步長(zhǎng)為0.2，, x y取值范圍-10到+10，圖(3-3)中Z比例為10，, x y都以步長(zhǎng)0.05變化，, x y取值范圍-5到+5。 -10 -8-6-4-20246810-10-8-6-4-20246810(a)(b)空間分布等值線(xiàn)分布圖 3-2 式(3-122)值的變化及等值線(xiàn)的分布 -5-4-3-2-1012345-5-4-3-2-1012345(a)(b)空間分布等值線(xiàn)分布圖 3-3 式(3-123)值的變化及等值線(xiàn)的分布 結(jié)果表示具有波動(dòng)性，表明床面有可能出現(xiàn)沙紋或沙波的變化形式，同時(shí)床面的變化形式與揚(yáng)沙含沙量的變化相關(guān)。圖中表示的床

36、面比較規(guī)則，但對(duì)復(fù)雜的水動(dòng)力條件和床面初始形態(tài)，床面變化具有擬周期性，同時(shí)變化程度非常復(fù)雜不規(guī)則。由于泥沙問(wèn)題的復(fù)雜性，目前對(duì)泥沙運(yùn)動(dòng)機(jī)理還沒(méi)有完全認(rèn)識(shí)清楚?，F(xiàn)有的各家對(duì)基本問(wèn)題的認(rèn)識(shí)及處理方法也不盡相同，甚至相差甚遠(yuǎn)。在20世紀(jì)70年代開(kāi)始出現(xiàn)了大量根據(jù)實(shí)測(cè)資料，利用計(jì)算機(jī)進(jìn)行回歸方法計(jì)算得出的阻力和輸沙力經(jīng)驗(yàn)公式，由于有大量的實(shí)驗(yàn)及天然資料做基礎(chǔ)，各種變量的覆蓋范圍大，這樣得出的一些公式精度很高，甚至超出了許多理論公式，其原因在于現(xiàn)有的理論體系往往不能夠把所有的泥沙的輸移過(guò)程完整地、正確地反映出來(lái)。本節(jié)內(nèi)容從泥沙擴(kuò)散的物理角度出發(fā)，基于現(xiàn)有的泥沙理論推導(dǎo)出兩個(gè)可分別計(jì)算輸沙和揚(yáng)沙的公式，在

37、理論擴(kuò)展、公式推導(dǎo)上做了一新的改進(jìn)，從公式結(jié)構(gòu)上看與現(xiàn)有懸移質(zhì)泥沙連續(xù)方程結(jié)構(gòu)相同，如將輸沙連續(xù)方程與揚(yáng)沙連續(xù)方程相疊加得到的全沙連續(xù)方程與現(xiàn)有懸移質(zhì)泥沙連續(xù)方程具有完全相同的形式，說(shuō)明方程推導(dǎo)正確，只是在劃分輸移泥沙和揚(yáng)動(dòng)泥沙標(biāo)準(zhǔn)上可能還存在進(jìn)一步研究的空間。3.2 初值條件及邊界條件的處理 初值條件可以由初始時(shí)刻的實(shí)測(cè)資料給出， 但所需條件在實(shí)際計(jì)算中一般難以全部獲得，只能通過(guò)估算加以補(bǔ)足，當(dāng)然估算得越接近實(shí)際越好。 對(duì)計(jì)算中所用的初始水流條件，一般以下游水位為起點(diǎn)，根據(jù)平均河床比降給出各點(diǎn)水位的初值，并令初始流速u(mài)、v均為0。含沙量初值有實(shí)測(cè)值或取為042。 通常對(duì)水流而言，進(jìn)口邊界給出

38、的流量過(guò)程線(xiàn)，出口邊界給定水位流量過(guò)程線(xiàn)。對(duì)懸移質(zhì)而言，進(jìn)口一般給定含沙量過(guò)程及相應(yīng)泥沙級(jí)配，出口認(rèn)為/0kSx 。 對(duì)于岸邊界條件，一般采用水流無(wú)滑移條件，即取u、v均為0。對(duì)懸移質(zhì)而言，岸邊界處取/0kSn。(n為邊界法線(xiàn)方向) 上游進(jìn)口邊界： 000( )( )( )uu tvv tSS t (3-124) 下游出口邊界： 00( )( )/0kuu tvv tSx (3-125) 3.3網(wǎng)格劃分網(wǎng)格的確定對(duì)模擬實(shí)際流動(dòng)的可靠性同樣重要。建網(wǎng)基本要求是：符合流動(dòng)特點(diǎn)；易于建立；比較光滑和規(guī)則；滿(mǎn)足精度和計(jì)算的穩(wěn)定要求；便于組成節(jié)約、高效的數(shù)據(jù)結(jié)構(gòu)，必要時(shí)可隨時(shí)依解的梯度作適應(yīng)性調(diào)整。本文

39、為了更好的模擬計(jì)算域內(nèi)的實(shí)際地形、地貌，采用的是無(wú)結(jié)構(gòu)有限體積網(wǎng)格劃分計(jì)算區(qū)域，使模型地形、地物與實(shí)際情況基本一致。一般而言，無(wú)結(jié)構(gòu)網(wǎng)格劃分采用三角形網(wǎng)格，因其有較為成熟的自動(dòng)生成技術(shù)，故自動(dòng)化程度較高，但是針對(duì)本文設(shè)計(jì)的二維河道情況，采用的是不規(guī)則四邊形組成的網(wǎng)格，主要考慮到：一方面若完全采用三角形網(wǎng)格作均勻劃分，反而會(huì)難以客觀(guān)反映小水量時(shí)水流沿主槽流動(dòng)的特點(diǎn)，有可能造成計(jì)算結(jié)果失真；另一方面，網(wǎng)格呈四邊形格式，局部輔以三角形網(wǎng)格，這樣不但可以很好地?cái)M合邊界形狀，而且格子數(shù)及其邊數(shù)，要比完全采用三角形網(wǎng)格時(shí)少，同時(shí)實(shí)際空間步長(zhǎng)也可以取大些。按照有限體積法布置二維網(wǎng)格的方式如圖 3-4 所示。

40、它取單元網(wǎng)格為控制體，在網(wǎng)格中心處計(jì)算水位H，在網(wǎng)格周邊通道的中點(diǎn)處計(jì)算流量Q。即在平衡計(jì)算時(shí)，沿控制體每一邊的法向通量用該邊中點(diǎn)處的通量作代表，乘以邊長(zhǎng)即為通量沿該邊的積分。中點(diǎn)的通量可用中心格式(如取相鄰兩格子形心處通量的平均)或逆風(fēng)格式確定。 計(jì)算所采用的方程改寫(xiě)成矢量形式， 按照有限體積法，將其在控制體內(nèi)進(jìn)行積分，對(duì)水位H和流量Q按時(shí)間交錯(cuò)方式計(jì)算(如圖 3-5)。 節(jié) 點(diǎn)通 道網(wǎng) 格通 道網(wǎng) 格T-2dtT+2dtT+dtT-dtT 圖 3-4 H和Q的空間布置方式 圖 3-5 H和Q的時(shí)間交錯(cuò)計(jì)算方式 3.4 水流挾沙能力1.懸移質(zhì)輸沙率公式的選取本文中采用的積分求輸沙率的方法，選

41、用了愛(ài)因斯坦挾沙力理論。*1230.211.62.303lg()sSvashgu aSIIk (3-131) 其中： 11110.216(1)ZZZAAIdA (3-132) 11210.216 ln(1)ZZZAAIdA (3-133) 1I和2I均是A及Z的積分函數(shù)， 可通過(guò)數(shù)值積分求解預(yù)先做成圖標(biāo)， 如圖3-7所示。 如果沒(méi)有床面形態(tài)時(shí)， 上述公式中的沙粒阻力對(duì)應(yīng)的剪切流速*u就等于*u。 愛(ài)因斯坦假定推移質(zhì)的運(yùn)動(dòng)速度為*11.6u，推移質(zhì)厚度為二倍粒徑。若求解懸移質(zhì)輸沙率時(shí)把沿垂線(xiàn)積分的下限取為2aD，記該處泥沙的體積比參考濃度vaS，則 */(11.62)vabsSguD (3-13

42、4) 2.泥沙計(jì)算輔助公式 (1)Manning公式48 21321uR Jn (3-135) 式中：n為Manning系數(shù)，也稱(chēng)為糙率，R是水力半徑，J是水力坡降。 (2)Manning-Strickler公式 1 66524()uRDggR J (3-136) 式中：u是斷面平均流速，R是沙粒阻力對(duì)應(yīng)的水力半徑，J同上，65D是泥沙級(jí)配中65%對(duì)應(yīng)的泥沙粒徑。 (3)泥沙顆粒統(tǒng)一沉速公式 29(9)sgDDD (3-137) 式中：v是水流的運(yùn)動(dòng)粘性，D是泥沙顆粒粒徑。 (7)起動(dòng)臨界流速 參考沙莫夫公式49 1/61.47()chUgDD (3-141) (8)止動(dòng)臨界流速 參考沙玉清公

43、式 3/40.011 0.39sUD (3-142) 1 21 21()()bssggD 式中：為推移質(zhì)強(qiáng)度。 第四章第四章 平面二維水流泥沙數(shù)學(xué)模型的數(shù)值模擬平面二維水流泥沙數(shù)學(xué)模型的數(shù)值模擬本文采用有限體積法離散方程， 它將計(jì)算區(qū)域劃分為網(wǎng)格，并使每個(gè)網(wǎng)格點(diǎn)周?chē)幸粋€(gè)互不重復(fù)的控制體；將待解微分方程對(duì)每一個(gè)控制體積分，從而得到一組離散方程，求解差分方程組，其求解就作為微分方程定解問(wèn)題的近似解。 4.1 水流控制方程的離散2122kTdtTTdtiiikikinidthhQ LdtqA (4-9) 式中，iA為第i個(gè)網(wǎng)格單元面積，ikL為i號(hào)網(wǎng)格的第K號(hào)通道長(zhǎng)度；ikQ為i號(hào)網(wǎng)格的第k號(hào)通道

44、的單寬流量。 動(dòng)量方程的求解，采用分類(lèi)簡(jiǎn)化處理的方法，把它概括為地面型通道、河網(wǎng)型通道和特殊通道來(lái)分類(lèi)處理。假定同一時(shí)段內(nèi)同一網(wǎng)格水位變化不大，同時(shí)由于平面淺水中對(duì)流項(xiàng)作用較小，在計(jì)算中略去不計(jì)，此外，有根據(jù)堰流流量公式計(jì)算的連續(xù)堤或缺口通道。1.河道型通道，即通道兩側(cè)網(wǎng)格均為河道型網(wǎng)格，動(dòng)量方程中保留局地加速度項(xiàng)、重力項(xiàng)和阻力項(xiàng)，利用差分方法離散得到河道型通道的動(dòng)量離散方程 2.地面型通道，即通道兩側(cè)單元為陸地地面，且通道上沒(méi)有堤防等阻水建筑物。考慮到庫(kù)區(qū)內(nèi)的地形變化，地面洪水演進(jìn)主要受到重力和阻力的作用，利用差分方法離散得到地面型通道的動(dòng)量離散方程 3.對(duì)于高于地面的阻水建筑物，如堤防、鐵

45、路、公路等，可以概化成連續(xù)堤或缺口堤通道，其流量采用寬頂堰溢流公式來(lái)計(jì)算，離散后得到 泥沙連續(xù)方程離散方法與水流連續(xù)方程離散方法相似，此處不再贅述。4.2 泥沙控制方程的離散 泥沙連續(xù)方程離散方法與水流連續(xù)方程離散方法相似，此處不再贅述。4.3 河床縱向變形方程的離散 對(duì)式子(3-104)的離散如下 22*TdtTbibidtSZZ (4-32) 式中：bZ是河底高程。 在計(jì)算河床變形要注意以下情況 SSSSSSSSSSStZb|, 0 |, 0 , 0 0, 0 且如果且如果但不滿(mǎn)足泥沙起動(dòng)條件如果且滿(mǎn)足泥沙起動(dòng)條件如果 (4-33) 4.4 模型求解過(guò)程 模型求解的過(guò)程具體步驟如下： (1

46、)給定計(jì)算所需要的基本參數(shù)，比如, ,xyn 等，代入初始條件和邊界條件； (2)對(duì)整個(gè)計(jì)算區(qū)域賦初值； (3)在時(shí)間步長(zhǎng)之內(nèi)，求解速度, u v，水深h； (4)根據(jù) Einstein 懸移質(zhì)理論計(jì)算vaS；應(yīng)用輸沙含沙量方程、揚(yáng)沙含沙量分別計(jì)算,S S、此后應(yīng)用懸浮指標(biāo)z及王尚毅的懸浮功理論判斷懸沙、底沙的交換比例，最后計(jì)算河床縱向變形bZ。 (5)至此一個(gè)時(shí)間步長(zhǎng)求解完畢，所求得的結(jié)果作為新的初始值，進(jìn)入下一個(gè)步長(zhǎng)的計(jì)算過(guò)程，重復(fù)第(3)(5)步驟，直至?xí)r間步長(zhǎng)循環(huán)到最終結(jié)果。 (6)輸出所需的計(jì)算結(jié)果。 第五章第五章 水流泥沙數(shù)學(xué)模型驗(yàn)證及應(yīng)用水流泥沙數(shù)學(xué)模型驗(yàn)證及應(yīng)用本模型的建立主要

47、用于廟宮水利樞紐庫(kù)區(qū)泥沙變化的研究，選擇了庫(kù)區(qū)河道為模擬對(duì)象，在分析庫(kù)區(qū)河道水流泥沙運(yùn)動(dòng)特點(diǎn)的基礎(chǔ)上，建立了二維水流泥沙數(shù)學(xué)模型，并用實(shí)測(cè)水文資料對(duì)模型進(jìn)行驗(yàn)證分析，在此基礎(chǔ)之上，進(jìn)行了水沙組合的計(jì)算。圖1-1廟宮水庫(kù)位置示意圖灤 平張 百 彎波 羅 諾風(fēng) 山灤河白 虎 溝成 子燕 上道 壩 子大 喚 起棋 盤(pán) 山龍 頭 山圍 場(chǎng)蟻螞吐河隆 化承 德廟 宮 水 庫(kù)釣 魚(yú) 臺(tái) 水 庫(kù)武烈河三 家荒 地大 廟伊 遜 河 流 域 平 面 位 置 圖廟宮水庫(kù)是河北省大型水庫(kù)中淤積浸沒(méi)最嚴(yán)重、問(wèn)題最突出、社會(huì)影響最大的水庫(kù)， 在海河流域乃至全國(guó)大中型水庫(kù)中也是淤積浸沒(méi)最嚴(yán)重的水庫(kù)之一。該水庫(kù)位于灤河的多

48、沙支流伊遜河上，由于建庫(kù)設(shè)計(jì)時(shí)資料短缺，上游來(lái)沙量計(jì)算偏小，對(duì)泥沙淤積的嚴(yán)重性估計(jì)不足，因此在水庫(kù)泄水建筑物的布置和運(yùn)用方式的擬定上，沒(méi)有充分考慮排沙的要求，致使大量泥沙淤積在庫(kù)區(qū)。水庫(kù)淤積形態(tài)為錐體淤積，壩前淤積高程 768.5m，超過(guò)死水位 0.5m。庫(kù)區(qū)河道比降由建庫(kù)前的 3.74變緩為 1.8，淤積末端距大壩 13.7km，上延系數(shù) 1.55。隨著淤積末端不斷上延，上游河床逐年抬高，距大壩 9.6km 的四合永公路橋處河床淤高 3m 多，使河床由原來(lái)的地面以下到現(xiàn)在與兩岸持平，甚至高出兩岸耕地及村基，產(chǎn)生嚴(yán)重浸沒(méi)。 5.1模型基本特征年際間輸沙量變化懸殊，廟宮水庫(kù)入庫(kù)最大沙量1360萬(wàn)

49、噸(1959年)，最小沙量61萬(wàn)噸(1981年)，兩者相差22.3倍。蟻螞吐河下河南站1959年輸沙量為1080萬(wàn)噸，1966年輸沙量?jī)H32萬(wàn)噸，變幅達(dá)33.8倍。水沙特點(diǎn)表現(xiàn)為水沙過(guò)程的不完全對(duì)應(yīng)一致，有時(shí)沙峰先于洪峰，有時(shí)反之；有時(shí)大水大沙、小水小沙、也有時(shí)大水小沙、小水大沙，主要是由于暴雨強(qiáng)度、過(guò)程和中心位置的差異以及各支流水土流失程度不同等原因造成的。5.2網(wǎng)格劃分廟宮水庫(kù)數(shù)學(xué)模型自水庫(kù)大壩向上游延伸，在四合永鎮(zhèn)北河道岔口分為兩支，一支沿伊遜河向西北至掌字水文站，另一支沿不澄河向東北至邊墻山水文站。模型主體南北走向，全長(zhǎng)約12.6km，模型最寬河段約1km。模型內(nèi)包含27個(gè)實(shí)測(cè)斷面，公

50、路橋、鐵路橋各一座。模型網(wǎng)格劃分以實(shí)測(cè)斷面為基礎(chǔ)，采用無(wú)結(jié)構(gòu)有限體積方法劃分網(wǎng)格，考慮主河槽走向，劃分后得到斷面194個(gè)，單元3451個(gè)，結(jié)點(diǎn)3626個(gè)，通道7076個(gè)。斷面平均間距為65m，斷面內(nèi)空間步長(zhǎng)為50m，單元平均面積為2124m2，模型覆蓋水域面積約為8km2，模型內(nèi)實(shí)測(cè)斷面分布見(jiàn)圖5-3，模型劃分網(wǎng)格見(jiàn)圖5-4，網(wǎng)格單元編號(hào)、通道編號(hào)、節(jié)點(diǎn)號(hào)的分布形式見(jiàn)圖5-5。 1234567891011121314151617181920212223242526272829303132332345678910111213141516171819202122232425262728293031

51、32333435363738123456789101112131415單 元 號(hào)通 道 號(hào)節(jié) 點(diǎn) 號(hào)：紅 色黑 色紫 色 5.3模型驗(yàn)證 1.水流模型驗(yàn)證 驗(yàn)證水流模型上游邊界采用伊遜河掌字水文站和不澄河邊墻山水文站1990年7月14日至15日歷時(shí)6小時(shí)流量過(guò)程（見(jiàn)圖5-5），下游邊界采用廟宮水庫(kù)水位流量關(guān)系過(guò)程（見(jiàn)表5-1），驗(yàn)證計(jì)算中沒(méi)有考慮排沙洞泄洪影響。假定大沽基面與國(guó)家高程基準(zhǔn)面之間的換算關(guān)系為：假定大沽基面-16.8681985國(guó)家高程基準(zhǔn)。水流數(shù)學(xué)模型驗(yàn)證伊遜河掌水文字站和不澄河邊墻山水文站1990年7月14日至15日歷時(shí)6小時(shí)水位過(guò)程，驗(yàn)證結(jié)果見(jiàn)圖5-8，驗(yàn)證計(jì)算流場(chǎng)見(jiàn)圖5-9

52、。從水位驗(yàn)證過(guò)程來(lái)看，計(jì)算值與實(shí)測(cè)值變化趨勢(shì)一致，數(shù)值偏差不大，基本能保證模型預(yù)報(bào)計(jì)算精度。從不同時(shí)刻流場(chǎng)上看，當(dāng)上游邊界產(chǎn)生大流量時(shí)，公路橋附近收縮斷面出現(xiàn)較大流速，下游大壩附近還沒(méi)有達(dá)到溢流水位；當(dāng)上游邊界流量回落時(shí)，公路橋附近斷面流速減弱，上游流量已傳遞到大壩附近，達(dá)到溢流水位，在壩前和庫(kù)區(qū)彎道附近產(chǎn)生較大流速，基本符合水流運(yùn)動(dòng)規(guī)律，說(shuō)明水流計(jì)算結(jié)果合理，水流模型可用于模型預(yù)報(bào)計(jì)算。從水位驗(yàn)證過(guò)程來(lái)看，計(jì)算值與實(shí)測(cè)值變化趨勢(shì)一致，數(shù)值偏差不大，基本能保證模型預(yù)報(bào)計(jì)算精度。從不同時(shí)刻流場(chǎng)上看，當(dāng)上游邊界產(chǎn)生大流量時(shí)，公路橋附近收縮斷面出現(xiàn)較大流速，下游大壩附近還沒(méi)有達(dá)到溢流水位；當(dāng)上游邊界

53、流量回落時(shí)，公路橋附近斷面流速減弱，上游流量已傳遞到大壩附近，達(dá)到溢流水位，在壩前和庫(kù)區(qū)彎道附近產(chǎn)生較大流速，基本符合水流運(yùn)動(dòng)規(guī)律，說(shuō)明水流計(jì)算結(jié)果合理，水流模型可用于模型預(yù)報(bào)計(jì)算2.泥沙模型驗(yàn)證泥沙數(shù)學(xué)模型驗(yàn)證采用計(jì)算模式驗(yàn)證的方法進(jìn)行，驗(yàn)證1990年汛期前后平均灘面的變化，汛期上游來(lái)流過(guò)程見(jiàn)圖5-6，上游邊界含沙量過(guò)程見(jiàn)圖5-10驗(yàn)證結(jié)果見(jiàn)圖5-11。20:0021:0022:0023:000:001:002:0090年7月14日7月15日時(shí)間（h）含沙量（kg/m ）305010015020020004000600080001000077077578078590年汛前平均床面90年汛后實(shí)

54、測(cè)平均床面90年汛后計(jì)算平均床面距壩里程（m）高程（m）34 567891011 12 13 14 1516 17 18 19 20 212223實(shí)測(cè)斷面序號(hào)圖 5-10 汛期上游邊界含沙量過(guò)程圖5-11 1990年汛期前后水庫(kù)平均灘面變化5.4 模型預(yù)報(bào)近年來(lái)隨著水環(huán)境關(guān)注意識(shí)的提高和水資源高效利用思路的確立，如何在水庫(kù)保水和水庫(kù)排沙間選擇優(yōu)化的水庫(kù)運(yùn)行模式是需要進(jìn)行研究和解決的重要問(wèn)題。根據(jù)廟宮水庫(kù)設(shè)計(jì)正常蓄水位778m(大沽高程，以下不加說(shuō)明處均采用大沽高程)和庫(kù)區(qū)死水位768m的要求，優(yōu)選即能泄水排沙又能保水興利的調(diào)水調(diào)沙方案是本次應(yīng)用的主要目的。研究過(guò)程擬采用庫(kù)區(qū)限蓄水位768m、7

55、70m、772m、774m、776m、778m六級(jí)水位標(biāo)準(zhǔn)，來(lái)水來(lái)沙過(guò)程采用19581959、19711975、19851989和1990的四個(gè)系列三種組合來(lái)模擬不同年份水沙過(guò)程，在18個(gè)計(jì)算系列中進(jìn)行方案優(yōu)選。下游排沙控制水位分別選取768m、770m、772m、774m、776m和778m，與三種水文類(lèi)型年組合系列再組合形成18種計(jì)算方案，方案命名以水文類(lèi)型年組合序號(hào)為先導(dǎo)字符，用橫線(xiàn)連接排沙控制水位值，計(jì)算方案見(jiàn)表6-2。5.5排沙量和淤積量分析根據(jù)方案的計(jì)算結(jié)果，得到沿廟宮水庫(kù)上游至下游各過(guò)流計(jì)算斷面的床面平均變形量曲線(xiàn)，便于比較將同一水文系列的情況在同一圖中給出(圖5-33至圖5-3

56、5)。從圖5-33至圖5-35中可以看出，三種不同類(lèi)型的來(lái)水來(lái)沙條件所塑造的河床變化形態(tài)相差較小，變化趨勢(shì)完全一致，只是下游不同的排沙控制水位對(duì)河床沖淤趨勢(shì)和沖淤量影響較大，故只給出水文條件下的最低排沙控制水位768以及最高排沙控制水位778情況下的庫(kù)區(qū)含沙量、沖淤變化分布、流量沿程變化(見(jiàn)圖5-36至圖5-45)。其他情況與此類(lèi)似，應(yīng)當(dāng)注意的是，輸沙含沙量與挾沙能力都是正值，而揚(yáng)沙含沙量、沖淤量都存在正負(fù)，對(duì)于揚(yáng)沙含沙量正值代表沖刷，負(fù)值代表淤積；沖淤量正值代表淤積，負(fù)值代表沖刷。從圖5-33至5-35可以看出，不同的排沙控制水位，對(duì)于離壩距離大于5200m(鐵路橋附近)的地方，影響很小，而影響范圍，主要集中在大壩到鐵路橋這一區(qū)段，變化趨勢(shì)是，隨著控制水位的不斷增大，排沙洞到大壩這一段淤積增加；排沙洞附近沖刷能力加強(qiáng)；排沙洞上游淤積增加?？傊褪?，排沙控制水位高，離排沙洞較遠(yuǎn)的地方淤積就會(huì)加重，而排沙洞附近的沖刷反而會(huì)提高。進(jìn)一步比較各方案之間在整體河道上的排沙效果，在表5-6中給出了廟宮水庫(kù)在三種不同類(lèi)型的來(lái)水來(lái)沙條件和不同排沙控制水位下各計(jì)算方案比較，13年水沙組合的來(lái)水量約為21470