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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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