摘要：采用有限元方法建立起一套河道平面二維水流泥 沙數學模型。在前人研究的基礎上，采用了質量集中的處理方法，提出了壓縮存儲的方法， 從而大大減少了計算存儲量。針對有限元法時間步長需取得較短問題，采用了“預報-校正-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解決工程實際問題，又大大減少了計算量。作者以下荊江監利河段為例進行泥沙沖淤計算，計算結果與實測值符合較好，從而證明了模型的可靠性。

關鍵詞：水流泥沙 有限元 模型驗證

三峽工程建成后，水庫將攔蓄大量泥沙，下泄水流含沙量減小，對三峽工程壩下游河道將產生以沖刷為主的影響，包括對荊江河段的河勢及荊江大堤帶來影響。為研究壩下游重點河段的河床沖淤分布、河勢變化、近岸流速變化等問題，一維模型顯得無能為力，但可采用平面二維模型來解決。有限元方法可采用無結構化網格，能很好地模擬不規則的幾何形狀，因此很適合于對天然河道的模擬。

然而，正如其它方法一樣，有限元法也有它的缺點，主要是計算存儲量和運算量較大。為揚長避短，使有限元方法能運用到對天然河道的模擬上來，本模型運用質量集中［4］的方法將系數矩陣轉化為三對角矩陣，并提出了緊湊的分塊壓縮存儲方法，從而大大減少了計算存儲量，使得計算能在一般微機上進行。采用質量集中方法的不足之處是時間步長需取得較短，且在河道模擬中尤為突出(因河道比較窄長，網格需劃分很細，而該法的穩定性要 求時間步長與網格尺度成正比)。針對該問題，筆者采用了“預報-校正-迭代［5］”的算法，該法可加大時間步長，同時有效避免了數值震蕩。針對長系列水沙條件下計算量較大問題，作者又提出了“非恒定-恒定-非恒定流”的算法，該算法既能解決工程實際問 題，又大大減少了計算量，使有限元方法能夠很好地運用于河道水流泥沙問題的實際計算。

1 基本方程

平面二維水流方程

(1)

(2)

(3)

懸移質泥沙擴散方程 

(4)

推移質不平衡輸移方程［6］

(5)

河床變形方程

由懸移質引起的河床變形方程為

(6)

由推移質引起的河床變形方程為

(7)

以上各式中U，V分別為垂線平均流速在x，y方向上的分量；Zs、Zb和H分別為水位、河底高程和水深；g為重力加速度；vt為水流紊動粘性系數；ρ為水的密度；τx、τy、舄瓂分別為底部切應力在x和y方(τx、τy)=，向上的分量：C為謝才系數，常用曼寧公式計算：C=H1/6/n；S和S*分別為垂線平均含沙量和挾沙力；N和N*分別為推移質輸沙量和推移質輸沙能力折算成全水深的泥沙濃度；εs為泥沙紊動擴散系數；ω為泥沙沉速；γ′為床沙干容重；α為懸移質泥沙恢復飽和系數，淤積時取0.25，沖刷時取1.0；β為推移質泥沙恢復飽和系數，取0.25。

懸移質泥沙共分成8組，水流挾沙力和分組挾沙力級配采用李義天方法［7］進行計算。

2有限元方程

將整個計算域剖分成一個三角形網格系統。每個三角形為一個單元，其編號為e，e=1，2，3， …，NE，NE為單元總數。單元三個頂點為節點，其局部編號為j=1，2，3(以逆時針為序)。節點的整體編號為i，i=1，2，3，…，NP，NP為節點總數。節點的整體編號與局部編號計算前一定要規定好。

引入插值函數：f=fi(t)φi(x，y)，φ為形函數。對方程(1)～(5)中的變量用插值函數近似表示，并使用伽遼金法［1］對方程進行整理變形，可得到積分方程

(8)

(9)

（10）

(11)

(12)

由此，經整理得到如下有限元方程

(13)

(14)

(15)

(16)

(17)

其中

i，j，k=1，2，3，…，NP.

3數值解法

為書寫方便，采用通用變量P來代替方程(13)～(17)中的變量Zs、U、V、S和N，用FP表示方程中等號右邊項，并用對角矩陣 來代替Aij，則方程可化成統一的形式： dPi/dt=FPi。

對該方程的求解，模型采用了“預報-校正-迭代［5］”的計算方法。其方法為：用二階顯式Adams公式作為預測公式，梯形公式(隱式)作校正公式，可構造Adams二階PC公式。離散后的方程不需聯接，計算過程穩定性好，時間步長可取得較長，同時還有效避免了數 值振蕩。具體計算過程如下

預報(18)

校正(19)

迭代給定誤差洌雜謁械1≤i≤NP，若，則令： ，否則令：= ；轉到校正，繼續迭代，直到滿足精度要求為止。當求出Hi、Ui、Vi、Si和Ni后，代入河床變形方程便可求得沖淤變形后的河床高程Zn+1bi，這樣便完成了一個時段的計算。以上角標n表時段，*表預報值，**表校正值，無*標記的表示該時段的計算結果。

在水位求解過程中還采用了Kawahara“選擇系數集中［4］”的經驗處理方法，引入了選擇性集中系數： α=0.8～1.0， 其作用相當于對計算結果進行“光滑”處理。

由于系數矩陣為NP×NP階矩陣，若直接進行存儲會占用大量的內存，從而導致計算無法進行。筆者發現系數矩陣中絕大部分為零元素，為此提出了壓縮存儲的方法，即只對系數矩陣中的非零元素進行存儲和計算，同時用輔助數組指明相應非零元素所在的位置。

對于長系列水沙條件的計算，提出了“非恒定-定非-定流”的算法，即將整個非恒定流量過程概化為梯級恒定流，對每個梯級恒定流采用非恒定流方程，以時間步長為迭代參數進行迭代，直到得到恒定的流場。

3數模計算有關問題處理

4.1初始條件及邊界條件

初始條件：可由初始刻實測資料給出。計算所需的初始條件在實際計算中一般難以全部獲得，只能通過估算加以補足，當然估算得越接近實際越好。

邊界條件：進口給定流量、含沙量及其級配，出口給定水位；對于不滑動岸邊界，取U=0，V= 0；對于滑動岸邊界，取邊界法線方向流速分量為0。

對于活動邊界，采用了動邊界模擬技術：對每個膖時段，采用計算的水位及水深值判別和區分水域和陸域計算節點。對陸域計算節點，使其保持一較小富余水深(Hmin=0.001m)，并取其糙率為一個接近于無窮大(如1010)的正數。

4.2三角形網格劃分

基于網格生成的基本思想，提出了四邊形法。對于單一河道，首先根據河勢將河道剖分為若干大的四邊形，再分別將這些大的四邊形剖分為若干小的四邊形，最后將四邊形的較短對角線相連，即形成三角形網格。對于分汊河道，將每條汊道當作單一河道處理后再進行拼接，便可形成整個計算域的網格。根據該方法，我們編制了通用的計算程序，只需輸入少量的信息，計算機便能自動生成網格，并給出節點坐標和單元關聯信息，最后配以屏幕顯示及圖形自動繪制。

4.3紊動粘性系數

根據零方程紊流模型，紊動粘性系數由vtαU*H公式確定。其中U*為摩阻流速；α為常數，經調試確定。對于泥沙擴散系數εs ，可近似取ε=vt。

4.3糙率

通過實測資料反求，并根據河道中不同的部位分塊調試糙率。

4.4床沙級配計算方法

在泥沙沖淤頻繁的河段，由于水流與泥沙的相互作用，使某組泥沙發生沖刷時，另一組泥沙可能發生淤積。因此，床沙級配的計算應能反映泥沙沖淤交替過程。本模型將計算河段內河床泥沙概化成三層，即表層(交換層)、次表層(擾動層)和深層。當表層某組泥沙發生沖刷時，處在次表層的該組泥沙將受到擾動，根據擾動強度的大小，確定該層泥沙進入表層參加交換的量。

5計算實例

河段概況：監利河彎段位于下荊江，上起塔市驛，下止于沙夾邊，長約20km。該河段平灘河寬約1400m，最寬處約3200m(烏龜洲)，最窄處約1000m，平灘水位下的平均水深約11m。河段內有一高程為30m(黃海)左右的烏龜洲(洲長約7km，寬約2km)將河道分成左右兩汊。本河段內有姚圻腦水文站，該站多年平均流量為11.370m3／s，多年平均含沙量為1.121kg/m3。河床主要由粉質粘土、砂粘土和細砂組成，河床表層床沙中值粒徑一般在0.00～0.22mm之間。

圖1 計算河段河勢及網格布置圖 Fig.1 The calculated river reach and the grids

計算區域的選取及網格劃分：選取塔市驛至沙夾邊長約20km的河段作為計算河段。采用三角形網格，將整個計算河段劃分為100個斷面，6732個三角形網格單元，共有3500個節點，見圖1。

地形條件：以1993年10月實測河道地形為計算起始地形，并用1996年10月實測地形進行驗證。

計算時段及水沙條件：選用1993年10月～1996年10月水沙資料，按流量共劃分為120個計算時段。進口按姚圻腦水文站流量和分組含沙量給出，出口水位根據沙夾邊與姚圻腦水位相關關系給出。

計算結果分析：

圖2為姚圻腦計算水位與實測水位對比圖，由圖中可見計算水位與實測水位吻合較好。

圖3為不同流量下流場示意圖，由圖中可見大水漫灘，小水歸槽，水流平順，流場合理。

圖4為1996年10月計算和實測的河道地形對比圖，由圖中可見通過計算得到的地形等高線與實測河道地形等高線的位置和范圍基本吻合。

圖2計算與實測水位對比圖 Fig.2Comparison between calculated and measured water level

圖3不同流量下流場示意圖 Fig.3Flow fields for different discharges

圖4計算與實測河道地形比較圖 Fig.4Comparison between calculated and measured channel

從驗證時段的沖淤總量來看，河段內實際淤積總量為2820萬m3，計算淤積總量為2852萬m3，相對誤差僅為1.1%。由此可見，本模型能較好反映本河段的河床變形情況。

6結論

(1)對于有限單元法在河道平面二維水流泥沙計算中的應用作了一些成功的探討，建立了一套行之有效的解法。有限元法具有能夠很好地模擬復雜邊界和河道地形等優點，但若不經過處理，就需要大量的存儲單元和很長的運算時間。為克服這些缺點，作者采用了質量集中的方法，并提出了壓縮存儲的方法，使得計算存儲量大為減少。另外，還采用了“預報-校正-迭代”的時間推進算法，并提出了“非恒定-恒定-非恒定流”的算法，大大減少了計算量。

(2)針對有限元前期單元剖分工作量巨大的問題，提出并運用四邊形法開發了三角形網格自動生成系統。該系統不僅能適用于單一河道，而且能適用于分汊和支流出匯、入匯等各種復雜情況。

(3)為檢驗模型的可靠性，文中以下荊江的監利河段為例進行了驗證計算。計算結果表明該模型能較好模擬本河段的水流泥沙運動及河床變形情況。

(4)由于天然河道的復雜性，目前關于泥沙運動理論還不夠成熟，且驗證資料尚不充分，因此數模計算精度有待作進一步的提高。