河道平面二維水沙數(shù)學(xué)模型的有限元法子

1、河道平面二維水沙數(shù)學(xué)模型的有限元法子摘要：采納 有限元法子建立起一套河道平面二維水流泥 沙數(shù)學(xué)模型。在前人鉆研 的根基 上，采納 了質(zhì)量集中的處理法子，提出了收縮 存儲(chǔ)的法子， 從而大大減少了盤算 存儲(chǔ)量。針對(duì)有限元法光陰 步長(zhǎng)需取得較短問題，采納 了“預(yù)報(bào)-校對(duì)-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解決工程實(shí)際問題，又大大減少了盤算 量。作者以下荊江監(jiān)利河段為例進(jìn)行泥沙沖淤盤算 ，盤算 效果 與實(shí)測(cè)值符合較好，從而證明了模型的可靠性。 要害 詞：水流泥沙 有限元 模型驗(yàn)證 三峽工程建成后，水庫(kù)將攔蓄大宗泥沙，下泄水流含沙量減小，對(duì)三峽工程壩下游河道將產(chǎn)生以沖洗 為主的影

2、響，包孕對(duì)荊江河段的河勢(shì)及荊江大堤帶來(lái)影響。為鉆研 壩下游重點(diǎn)河段的河床沖淤散播、河勢(shì)變更 、近岸流速變更 等問題，一維模型顯得無(wú)能為力，但可采納 平面二維模型來(lái)解決。有限元法子可采納 無(wú)結(jié)構(gòu)化網(wǎng)格，能很好地模仿 不規(guī)矩 的幾何形狀，因此很適宜于對(duì)天然河道的模仿 。然而，正如其它法子一樣，有限元法也有它的錯(cuò)誤，首要是盤算 存儲(chǔ)量和運(yùn)算量較大。為揚(yáng)長(zhǎng)避短，使有限元法子能運(yùn)用到對(duì)天然河道的模仿 上來(lái)，本模型運(yùn)用質(zhì)量集中4的法子將系數(shù)矩陣轉(zhuǎn)化為三對(duì)角矩陣，并提出了緊湊的分塊收縮 存儲(chǔ)法子，從而大大減少了盤算 存儲(chǔ)量，使得盤算 能在一般微機(jī)上進(jìn)行。采納 質(zhì)量集中法子的不足之處是光陰 步長(zhǎng)需取得較短，且

3、在河道模仿 中尤為突出(因河道對(duì)比 窄長(zhǎng)，網(wǎng)格需劃分很細(xì)，而該法的穩(wěn)固 性要 求光陰 步長(zhǎng)與網(wǎng)格標(biāo)準(zhǔn) 成正比)。針對(duì)該問題，筆者采納 了“預(yù)報(bào)-校對(duì)-迭代5”的算法，該法可加大光陰 步長(zhǎng)，同時(shí)有效避免了數(shù)值震動(dòng) 。針對(duì)長(zhǎng)系列水沙條件下盤算 量較大問題，作者又提出了“非恒定-恒定-非恒定流”的算法，該算法既能解決工程實(shí)際問 題，又大大減少了盤算 量，使有限元法子能夠很好地運(yùn)用于河道水流泥沙問題的實(shí)際盤算 。1 根基方程平面二維水流方程(1) (2)(3)懸移質(zhì)泥沙擴(kuò)散方程 (4)推移質(zhì)不平衡輸移方程6(5)河床變形方程由懸移質(zhì)引起的河床變形方程為(6)由推移質(zhì)引起的河床變形方程為(7)以上各式中

4、U，V分辨為垂線平均流速在x,y方向上的分量；Zs、Zb和H分辨為水位、河底高程和水深；g為重力加速度；vt為水流紊動(dòng)粘性系數(shù)；為水的密度；x、y、舄瓂分辨為底部切應(yīng)力在x和y方(x、y)=，向上的分量：C為謝才系數(shù)，常用曼寧公式盤算 ：C=H1/6/n；S和S*分辨為垂線平均含沙量和挾沙力；N和N*分辨為推移質(zhì)輸沙量和推移質(zhì)輸沙能力折算成全水深的泥沙濃度；s為泥沙紊動(dòng)擴(kuò)散系數(shù)；為泥沙沉速；為床沙干容重；為懸移質(zhì)泥沙恢復(fù)飽和系數(shù)，淤積時(shí)取0.25，沖洗 時(shí)取1.0；為推移質(zhì)泥沙恢復(fù)飽和系數(shù)，取0.25。懸移質(zhì)泥沙共分成8組，水流挾沙力和分組挾沙力級(jí)配采納 李義天法子7進(jìn)行盤算 。2有限元方程

5、將全部 盤算 域剖分成一個(gè)三角形網(wǎng)格系統(tǒng)。每個(gè)三角形為一個(gè)單元，其編號(hào)為e,e=1,2,3， ，NE,NE為單元總數(shù)。單元三個(gè)頂點(diǎn)為節(jié)點(diǎn)，其局部編號(hào)為j=1,2,3(以逆時(shí)針為序)。節(jié)點(diǎn)的整體編號(hào)為i,i=1,2,3,，NP，NP為節(jié)點(diǎn)總數(shù)。節(jié)點(diǎn)的整體編號(hào)與局部編號(hào)盤算 前必然要規(guī)定好。引入插值函數(shù)：f=fi(t)i(x,y),為形函數(shù)。對(duì)方程(1)(5)中的變量用插值函數(shù)近似表現(xiàn) ，并應(yīng)用 伽遼金法1對(duì)方程進(jìn)行收拾 變形，可得到積分方程(8)(9) （10）(11)(12)由此，經(jīng)收拾 得到如下有限元方程(13)(14)(15)(16)(17)其中i,j,k=1,2,3,，NP.3數(shù)值解法為

6、書寫方便，采納 通用變量P來(lái)代替方程(13)(17)中的變量Zs、U、V、S和N，用FP表現(xiàn) 方程中等號(hào)右邊項(xiàng)，并用對(duì)角矩陣 來(lái)代替Aij，則方程可化成統(tǒng)一的情勢(shì) ： dPi/dt=FPi。對(duì)該方程的求解，模型采納 了“預(yù)報(bào)-校對(duì)-迭代5”的盤算 法子。其法子為：用二階顯式Adams公式作為預(yù)測(cè)公式，梯形公式(隱式)作校對(duì)公式，可結(jié)構(gòu) Adams二階PC公式。離散后的方程不需聯(lián)接，盤算 歷程 穩(wěn)固 性好，光陰 步長(zhǎng)可取得較長(zhǎng)，同時(shí)還有效避免了數(shù) 值振蕩。具體盤算 歷程 如下預(yù)報(bào)(18)校對(duì)(19)迭代給定誤差洌雜謁械1iNP，若，則令：，否則令：= ；轉(zhuǎn)到校對(duì)，持續(xù) 迭代，直到滿足精度請(qǐng)求 為

7、止。當(dāng)求出Hi、Ui、Vi、Si和Ni后，代入河床變形方程便可求得沖淤變形后的河床高程Zn+1bi，這樣便完成了一個(gè)時(shí)段的盤算 。以上角標(biāo)n表時(shí)段，*表預(yù)報(bào)值，*表校對(duì)值，無(wú)*標(biāo)記的表現(xiàn) 該時(shí)段的盤算 效果 。在水位求解歷程 中還采納 了Kawahara“選擇系數(shù)集中4”的經(jīng)驗(yàn)處理法子，引入了選擇性集中系數(shù)： =0.81.0， 其作用相當(dāng)于對(duì)盤算 效果 進(jìn)行“光滑”處理。 由于系數(shù)矩陣為NPNP階矩陣，若直接進(jìn)行存儲(chǔ)會(huì)占用大宗的內(nèi)存，從而導(dǎo)致盤算 無(wú)法進(jìn)行。筆者創(chuàng)造系數(shù)矩陣中絕大部分為零元素，為此提出了收縮 存儲(chǔ)的法子，即只對(duì)系數(shù)矩陣中的非零元素進(jìn)行存儲(chǔ)和盤算 ，同時(shí)用贊助數(shù)組指明相應(yīng)非零元素

8、所在的地位 。對(duì)于長(zhǎng)系列水沙條件的盤算 ，提出了“非恒定-定非-定流”的算法，即將全部 非恒定流量歷程 概化為梯級(jí)恒定流，對(duì)每個(gè)梯級(jí)恒定流采納 非恒定流方程，以光陰 步長(zhǎng)為迭代參數(shù)進(jìn)行迭代，直到得到恒定的流場(chǎng)。3數(shù)模盤算 有關(guān)問題處理4.1初始條件及邊界條件初始條件：可由初始刻實(shí)測(cè)材料 給出。盤算 所需的初始條件在實(shí)際盤算 中一般難以整個(gè) 獲得，只能通過估算加以補(bǔ)足，當(dāng)然估算得越接近實(shí)際越好。邊界條件：進(jìn)口給定流量、含沙量及其級(jí)配，出口給定水位；對(duì)于不滑動(dòng)岸邊界，取U=0,V= 0；對(duì)于滑動(dòng)岸邊界，取邊界法線方向流速分量為0。對(duì)于運(yùn)動(dòng) 邊界，采納 了動(dòng)邊界模仿 技巧 ：對(duì)每個(gè)膖時(shí)段，采納 盤算

9、 的水位及水深值判別和區(qū)分水域和陸域盤算 節(jié)點(diǎn)。對(duì)陸域盤算 節(jié)點(diǎn)，使其維持 一較小富裕水深(Hmin=0.001m)，并取其糙率為一個(gè)接近于無(wú)窮大(如1010)的正數(shù)。4.2三角形網(wǎng)格劃分基于網(wǎng)格生成的根基思想，提出了四邊形法。對(duì)于單一河道，首先根據(jù)河勢(shì)將河道剖分為若干大的四邊形，再分辨將這些大的四邊形剖分為若干小的四邊形，最后將四邊形的較短對(duì)角線相連，即形成三角形網(wǎng)格。對(duì)于分汊河道，將每條汊道當(dāng)作單一河道處理后再進(jìn)行拼接，便可形成全部 盤算 域的網(wǎng)格。根據(jù)該法子，我們編制了通用的盤算 程序，只需輸入少量的信息，盤算 機(jī)便能主動(dòng) 生成網(wǎng)格，并給出節(jié)點(diǎn)坐標(biāo)和單元關(guān)聯(lián)信息，最后配以屏幕顯示及圖形主

10、動(dòng) 繪制。4.3紊動(dòng)粘性系數(shù)根據(jù)零方程紊流模型，紊動(dòng)粘性系數(shù)由vtU*H公式斷定 。其中U*為摩阻流速；為常數(shù)，經(jīng)調(diào)試斷定 。對(duì)于泥沙擴(kuò)散系數(shù)s ,可近似取=vt。4.3糙率通過實(shí)測(cè)材料 反求，并根據(jù)河道中不同的部位分塊調(diào)試糙率。4.4床沙級(jí)配盤算 法子在泥沙沖淤頻繁的河段，由于水流與泥沙的相互作用，使某組泥沙產(chǎn)生 沖洗 時(shí)，另一組泥沙可能產(chǎn)生 淤積。因此，床沙級(jí)配的盤算 應(yīng)能反響泥沙沖淤交替歷程 。本模型將盤算 河段內(nèi)河床泥沙概化成三層，即表層(交換層)、次表層(擾動(dòng)層)和深層。當(dāng)表層某組泥沙產(chǎn)生 沖洗 時(shí)，處在次表層的該組泥沙將受到擾動(dòng)，根據(jù)擾動(dòng)強(qiáng)度的大小，斷定 該層泥沙進(jìn)入表層參加交換的

11、量。5盤算 實(shí)例河段概況：監(jiān)利河彎段位于下荊江，上起塔市驛，下止于沙夾邊，長(zhǎng)約20km。該河段平灘河寬約1400m，最寬處約3200m(烏龜洲)，最窄處約1000m，平灘水位下的平均水深約11m。河段內(nèi)有一高程為30m(黃海)左右的烏龜洲(洲長(zhǎng)約7km，寬約2km)將河道分成左右兩汊。本河段內(nèi)有姚圻腦水文站，該站多年平均流量為11.370m3s，多年平均含沙量為1.121kg/m3。河床首要由粉質(zhì)粘土、砂粘土和細(xì)砂組成，河床表層床沙中值粒徑一般在0.000.22mm之間。圖1 盤算 河段河勢(shì)及網(wǎng)格安排圖Fig.1 The calculated river reach and the grids

12、盤算 區(qū)域的選取及網(wǎng)格劃分：選取塔市驛至沙夾邊長(zhǎng)約20km的河段作為盤算 河段。采納 三角形網(wǎng)格，將全部 盤算 河段劃分為100個(gè)斷面，6732個(gè)三角形網(wǎng)格單元，共有3500個(gè)節(jié)點(diǎn)，見圖1。地形條件：以1993年10月實(shí)測(cè)河道地形為盤算 起始地形，并用1996年10月實(shí)測(cè)地形進(jìn)行驗(yàn)證。盤算 時(shí)段及水沙條件：選用1993年10月1996年10月水沙材料 ，按流量共劃分為120個(gè)盤算 時(shí)段。進(jìn)口按姚圻腦水文站流量和分組含沙量給出，出口水位根據(jù)沙夾邊與姚圻腦水位相干 關(guān)系給出。盤算 效果 分析：圖2為姚圻腦盤算 水位與實(shí)測(cè)水位對(duì)照 圖，由圖中可見盤算 水位與實(shí)測(cè)水位吻合較好。圖3為不同流量下流場(chǎng)示意

13、圖，由圖中可見大水漫灘，小水歸槽，水流平順，流場(chǎng)合理。圖4為1996年10月盤算 和實(shí)測(cè)的河道地形對(duì)照 圖，由圖中可見通過盤算 得到的地形等高線與實(shí)測(cè)河道地形等高線的地位 和領(lǐng)域根基吻合。圖2盤算 與實(shí)測(cè)水位對(duì)照 圖Fig.2Comparison between calculated andmeasured water level圖3不同流量下流場(chǎng)示意圖Fig.3Flow fields for different discharges圖4盤算 與實(shí)測(cè)河道地形對(duì)比 圖Fig.4Comparison between calculated and measured channel從驗(yàn)證時(shí)段的沖淤總量

14、來(lái)看，河段內(nèi)實(shí)際淤積總量為2820萬(wàn)m3，盤算 淤積總量為2852萬(wàn)m3，相對(duì)誤差僅為1.1%。由此可見，本模型能較好反響本河段的河床變形情況。6結(jié)論(1)對(duì)于有限單元法在河道平面二維水流泥沙盤算 中的利用 作了一些成功的探討，建立了一套行之有效的解法。有限元法具有能夠很好地模仿 繁雜 邊界和河道地形等優(yōu)點(diǎn)，但若不經(jīng)過處理，就需要大宗的存儲(chǔ)單元和很長(zhǎng)的運(yùn)算光陰 。為克服這些錯(cuò)誤，作者采納 了質(zhì)量集中的法子，并提出了收縮 存儲(chǔ)的法子，使得盤算 存儲(chǔ)量大為減少。另外，還采納 了“預(yù)報(bào)-校對(duì)-迭代”的光陰 推動(dòng) 算法，并提出了“非恒定-恒定-非恒定流”的算法，大大減少了盤算 量。(2)針對(duì)有限元前期

15、單元剖分工作量龐大的問題，提出并運(yùn)用四邊形法開發(fā)了三角形網(wǎng)格主動(dòng) 生成系統(tǒng)。該系統(tǒng)不僅能實(shí)用 于單一河道，而且能實(shí)用 于分汊和支流出匯、入?yún)R等各種繁雜 情況。(3)為檢驗(yàn)?zāi)Ｐ偷目煽啃?，文中以下荊江的監(jiān)利河段為例進(jìn)行了驗(yàn)證盤算 。盤算 效果 表明該模型能較好模仿 本河段的水流泥沙運(yùn)動(dòng)及河床變形情況。(4)由于天然河道的繁雜 性，目前關(guān)于泥沙運(yùn)動(dòng)理論還不夠成熟，且驗(yàn)證材料 尚不充沛 ，因此數(shù)模盤算 精度有待作進(jìn)一步的進(jìn)步 。參考文獻(xiàn) 1周雪漪.盤算 水力學(xué)M.清華大學(xué)出版社，1995.9. 2謝鑒衡主編.河流模仿 M.水利電力出版社，1988.6. 3張瑞瑾，謝鑒衡等.河流泥沙動(dòng)力學(xué)M.水利電力出版社，1989.6. 4Kawahara,M.,Hirano,H.,Tsubota,K.,andInayuki,K.SelectiveLumpingFiniteElementMethodforShallowWaterFlow,InternationJournalforNumericalMethodinFluids,1980,vol.2:89-112. 5韋直林.河道水流泥沙問題的一種有限元解法J.武漢水利電