水利工程論文-河道平面二維水沙數(shù)學(xué)模型的有限元方法摘要：采用有限元方法建立起一套河道平面二維水流泥沙數(shù)學(xué)模型。在前人研究的基礎(chǔ)上，采用了質(zhì)量集中的處理方法，提出了壓縮存儲(chǔ)的方法，從而大大減少了計(jì)算存儲(chǔ)量。針對(duì)有限元法時(shí)間步長(zhǎng)需取得較短問(wèn)題，采用了“預(yù)報(bào)-校正-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解決工程實(shí)際問(wèn)題，又大大減少了計(jì)算量。作者以下荊江監(jiān)利河段為例進(jìn)行泥沙沖淤計(jì)算，計(jì)算結(jié)果與實(shí)測(cè)值符合較好，從而證明了模型的可靠性。關(guān)鍵詞：水流泥沙有限元模型驗(yàn)證三峽工程建成后，水庫(kù)將攔蓄大量泥沙，下泄水流含沙量減小，對(duì)三峽工程壩下游河道將產(chǎn)生以沖刷為主的影響，包括對(duì)荊江河段的河勢(shì)及荊江大堤帶來(lái)影響。為研究壩下游重點(diǎn)河段的河床沖淤分布、河勢(shì)變化、近岸流速變化等問(wèn)題，一維模型顯得無(wú)能為力，但可采用平面二維模型來(lái)解決。有限元方法可采用無(wú)結(jié)構(gòu)化網(wǎng)格，能很好地模擬不規(guī)則的幾何形狀，因此很適合于對(duì)天然河道的模擬。然而，正如其它方法一樣，有限元法也有它的缺點(diǎn)，主要是計(jì)算存儲(chǔ)量和運(yùn)算量較大。為揚(yáng)長(zhǎng)避短，使有限元方法能運(yùn)用到對(duì)天然河道的模擬上來(lái)，本模型運(yùn)用質(zhì)量集中4的方法將系數(shù)矩陣轉(zhuǎn)化為三對(duì)角矩陣，并提出了緊湊的分塊壓縮存儲(chǔ)方法，從而大大減少了計(jì)算存儲(chǔ)量，使得計(jì)算能在一般微機(jī)上進(jìn)行。采用質(zhì)量集中方法的不足之處是時(shí)間步長(zhǎng)需取得較短，且在河道模擬中尤為突出(因河道比較窄長(zhǎng)，網(wǎng)格需劃分很細(xì)，而該法的穩(wěn)定性要求時(shí)間步長(zhǎng)與網(wǎng)格尺度成正比)。針對(duì)該問(wèn)題，筆者采用了“預(yù)報(bào)-校正-迭代5”的算法，該法可加大時(shí)間步長(zhǎng)，同時(shí)有效避免了數(shù)值震蕩。針對(duì)長(zhǎng)系列水沙條件下計(jì)算量較大問(wèn)題，作者又提出了“非恒定-恒定-非恒定流”的算法，該算法既能解決工程實(shí)際問(wèn)題，又大大減少了計(jì)算量，使有限元方法能夠很好地運(yùn)用于河道水流泥沙問(wèn)題的實(shí)際計(jì)算。1基本方程平面二維水流方程(1)(2)(3)(4)推移質(zhì)不平衡輸移方程6(5)河床變形方程由懸移質(zhì)引起的河床變形方程為(6)由推移質(zhì)引起的河床變形方程為(7)以上各式中U，V分別為垂線平均流速在x,y方向上的分量；Zs、Zb和H分別為水位、河底高程和水深；g為重力加速度；vt為水流紊動(dòng)粘性系數(shù)；為水的密度；x、y、舄瓂分別為底部切應(yīng)力在x和y方(x、y)=，向上的分量：C為謝才系數(shù)，常用曼寧公式計(jì)算：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)行計(jì)算。2有限元方程將整個(gè)計(jì)算域剖分成一個(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)計(jì)算前一定要規(guī)定好。引入插值函數(shù)：f=fi(t)i(x,y),為形函數(shù)。對(duì)方程(1)(5)中的變量用插值函數(shù)近似表示，并使用伽遼金法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ù)值解法為書寫方便，采用通用變量P來(lái)代替方程(13)(17)中的變量Zs、U、V、S和N，用FP表示方程中等號(hào)右邊項(xiàng)，并用對(duì)角矩陣來(lái)代替Aij，則方程可化成統(tǒng)一的形式：dPi/dt=FPi。對(duì)該方程的求解，模型采用了“預(yù)報(bào)-校正-迭代5”的計(jì)算方法。其方法為：用二階顯式Adams公式作為預(yù)測(cè)公式，梯形公式(隱式)作校正公式，可構(gòu)造Adams二階PC公式。離散后的方程不需聯(lián)接，計(jì)算過(guò)程穩(wěn)定性好，時(shí)間步長(zhǎng)可取得較長(zhǎng)，同時(shí)還有效避免了數(shù)值振蕩。具體計(jì)算過(guò)程如下預(yù)報(bào)(18)校正(19)迭代給定誤差洌雜謁械1iNP，若，則令：，否則令：=；轉(zhuǎn)到校正，繼續(xù)迭代，直到滿足精度要求為止。當(dāng)求出Hi、Ui、Vi、Si和Ni后，代入河床變形方程便可求得沖淤變形后的河床高程Zn+1bi，這樣便完成了一個(gè)時(shí)段的計(jì)算。以上角標(biāo)n表時(shí)段，*表預(yù)報(bào)值，*表校正值，無(wú)*標(biāo)記的表示該時(shí)段的計(jì)算結(jié)果。在水位求解過(guò)程中還采用了Kawahara“選擇系數(shù)集中4”的經(jīng)驗(yàn)處理方法，引入了選擇性集中系數(shù)：=0.81.0，其作用相當(dāng)于對(duì)計(jì)算結(jié)果進(jìn)行“光滑”處理。由于系數(shù)矩陣為NPNP階矩陣，若直接進(jìn)行存儲(chǔ)會(huì)占用大量的內(nèi)存，從而導(dǎo)致計(jì)算無(wú)法進(jìn)行。筆者發(fā)現(xiàn)系數(shù)矩陣中絕大部分為零元素，為此提出了壓縮存儲(chǔ)的方法，即只對(duì)系數(shù)矩陣中的非零元素進(jìn)行存儲(chǔ)和計(jì)算，同時(shí)用輔助數(shù)組指明相應(yīng)非零元素所在的位置。對(duì)于長(zhǎng)系列水沙條件的計(jì)算，提出了“非恒定-定非-定流”的算法，即將整個(gè)非恒定流量過(guò)程概化為梯級(jí)恒定流，對(duì)每個(gè)梯級(jí)恒定流采用非恒定流方程，以時(shí)間步長(zhǎng)為迭代參數(shù)進(jìn)行迭代，直到得到恒定的流場(chǎng)。3數(shù)模計(jì)算有關(guān)問(wèn)題處理4.1初始條件及邊界條件初始條件：可由初始刻實(shí)測(cè)資料給出。計(jì)算所需的初始條件在實(shí)際計(jì)算中一般難以全部獲得，只能通過(guò)估算加以補(bǔ)足，當(dāng)然估算得越接近實(shí)際越好。邊界條件：進(jìn)口給定流量、含沙量及其級(jí)配，出口給定水位；對(duì)于不滑動(dòng)岸邊界，取U=0,V=0；對(duì)于滑動(dòng)岸邊界，取邊界法線方向流速分量為0。對(duì)于活動(dòng)邊界，采用了動(dòng)邊界模擬技術(shù)：對(duì)每個(gè)膖時(shí)段，采用計(jì)算的水位及水深值判別和區(qū)分水域和陸域計(jì)算節(jié)點(diǎn)。對(duì)陸域計(jì)算節(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)行拼接，便可形成整個(gè)計(jì)算域的網(wǎng)格。根據(jù)該方法，我們編制了通用的計(jì)算程序，只需輸入少量的信息，計(jì)算機(jī)便能自動(dòng)生成網(wǎng)格，并給出節(jié)點(diǎn)坐標(biāo)和單元關(guān)聯(lián)信息，最后配以屏幕顯示及圖形自動(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糙率通過(guò)實(shí)測(cè)資料反求，并根據(jù)河道中不同的部位分塊調(diào)試糙率。4.4床沙級(jí)配計(jì)算方法在泥沙沖淤頻繁的河段，由于水流與泥沙的相互作用，使某組泥