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

1、河道平面二維水沙數(shù)學(xué)模型的有限元方法    然而，正如其它方法一樣，有限元法也有它的缺點(diǎn)，主要是計(jì)算存儲(chǔ)量和運(yùn)算量較大。為揚(yáng)長(zhǎng)避短，使有限元方法能運(yùn)用到對(duì)天然河道的模擬上來，本模型運(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ì)該問題，筆者采用了 “ 預(yù)報(bào) - 校正 - 迭代 5 ” 的算法，該法可加大時(shí)

2、間步長(zhǎng)，同時(shí)有效避免了數(shù)值震蕩。針對(duì)長(zhǎng)系列水沙條件下計(jì)算量較大問題，作者又提出了 “ 非恒定 - 恒定 - 非恒定流 ” 的算法，該算法既能解決工程實(shí)際問 題，又大大減少了計(jì)算量，使有限元方法能夠很好地運(yùn)用于河道水流泥沙問題的實(shí)際計(jì)算。 1 基本方程 平面二維水流方程 (1) (2) (3) 懸移質(zhì)泥沙擴(kuò)散方程 (4) 推移質(zhì)不平衡輸移方程6 (5) 河床變形方程 由懸移質(zhì)引起的河床變形方程為 <!endif> (6) 由推移質(zhì)引起的河床變形方程為 (7) 以上各式中 U ， V 分別為垂線平均流速在 x,y 方向上的分量； Zs 、 Zb 和 H 分別為水位、河底高程和水深； g

3、 為重力加速度； 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組，水流挾沙力和分組挾沙

4、力級(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)整理得到如下有限元方

5、程 (13) (14) (15) (16) (17) 其中 i,j,k=1,2,3, ， NP. 對(duì)該方程的求解，模型采用了 “ 預(yù)報(bào) - 校正 - 迭代 5 ” 的 計(jì)算 方法 。其方法為：用二階顯式 Adams 公式作為預(yù)測(cè)公式，梯形公式 ( 隱式 ) 作校正公式，可構(gòu)造 Adams 二階 PC 公式。離散后的方程不需聯(lián)接，計(jì)算過程穩(wěn)定性好，時(shí)間步長(zhǎng)可取得較長(zhǎng)，同時(shí)還有效避免了數(shù) 值振蕩。具體計(jì)算過程如下 預(yù)報(bào) (18) 校正 (19) 迭代給定誤差洌雜謁械 1iNP ，若 ，則令：      ，否則令： = ；轉(zhuǎn)到校正，繼續(xù)迭代，直到滿足

6、精度要求為止。當(dāng)求出 Hi 、 Ui 、 Vi 、 Si 和 Ni 后，代入 河床變形方程便可求得沖淤變形后的河床高程 Zn+1bi ，這樣便完成了一個(gè)時(shí)段的計(jì)算。以上角標(biāo) n 表時(shí)段， * 表預(yù)報(bào)值， * 表校正值，無 * 標(biāo)記的表示該時(shí)段的計(jì)算結(jié)果。     在水位求解過程中還采用了 Kawahara“ 選擇系數(shù)集中 4 ” 的 經(jīng)驗(yàn) 處理方法，引入了選擇性集中系數(shù)： =0.8 1.0 ，    其作用相當(dāng)于對(duì)計(jì)算結(jié)果進(jìn)行 “ 光滑 ” 處理。    由于系數(shù)矩陣為 NP

7、×NP 階矩陣，若直接進(jìn)行存儲(chǔ)會(huì)占用大量的內(nèi)存，從而導(dǎo)致計(jì)算無法進(jìn)行。筆者發(fā)現(xiàn)系數(shù)矩陣中絕大部分為零元素，為此提出了壓縮存儲(chǔ)的方法，即只對(duì)系數(shù)矩陣中的非零元素進(jìn)行存儲(chǔ)和計(jì)算，同時(shí)用輔助數(shù)組指明相應(yīng)非零元素所在的位置。 對(duì)于長(zhǎng)系列水沙條件的計(jì)算，提出了 “ 非恒定 - 定非 - 定流 ” 的算法，即將整個(gè)非恒定流量過程概化為梯級(jí)恒定流，對(duì)每個(gè)梯級(jí)恒定流采用非恒定流方程，以時(shí)間步長(zhǎng)為迭代參數(shù)進(jìn)行迭代，直到得到恒定的流場(chǎng)。 3 數(shù)模計(jì)算有關(guān) 問題 處理 4.1 初始條件及邊界條件 初始條件：可由初始刻實(shí)測(cè)資料給出。計(jì)算所需的初始條件在實(shí)際計(jì)算中一般難以全部獲得，只能通過估算加以補(bǔ)足，當(dāng)然估

8、算得越接近實(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è)接近于無窮大 ( 如 1010) 的正數(shù)。 4.2 三角形網(wǎng)格劃分 基于網(wǎng)格生成的基本思想，提出了四邊形法。對(duì)于單一河道，首先根據(jù)河勢(shì)將河道剖分為若干大的四邊形，再分別將這些大的四邊形剖分為若干小的四邊形，最后將四邊形的較短對(duì)角線相連，即形成三

9、角形網(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 糙率 通過實(shí)測(cè)資料反求，并根據(jù)河道中不同的部位分塊調(diào)試糙率。 4.4 床沙級(jí)配計(jì)算方法 在泥沙沖淤頻繁的河段，由于水流與泥沙的相互作用，使某組泥沙發(fā)生沖刷時(shí)，另一組泥沙可能發(fā)生淤

10、積。因此，床沙級(jí)配的計(jì)算應(yīng)能反映泥沙沖淤交替過程。本模型將計(jì)算河段內(nèi)河床泥沙概化成三層，即表層(交換層)、次表層(擾動(dòng)層)和深層。當(dāng)表層某組泥沙發(fā)生沖刷時(shí)，處在次表層的該組泥沙將受到擾動(dòng)，根據(jù)擾動(dòng)強(qiáng)度的大小，確定該層泥沙進(jìn)入表層參加交換的量。 5 計(jì)算實(shí)例 河段概況：監(jiān)利河彎段位于下荊江，上起塔市驛，下止于沙夾邊，長(zhǎng)約 20km 。該河段平灘河寬約 1400m ，最寬處約 3200m ( 烏龜洲 ) ，最窄處約 1000m ，平灘水位下的平均水深約 11m 。河段內(nèi)有一高程為 30m ( 黃海 ) 左右的烏龜洲 ( 洲長(zhǎng)約 7km ，寬約 2km ) 將河道分成左右兩汊。本河段內(nèi)有姚圻腦水文站

11、，該站多年平均流量為 11.370m3 s ，多年平均含沙量為 1.121kg /m3 。河床主要由粉質(zhì)粘土、砂粘土和細(xì)砂組成，河床表層床沙中值粒徑一般在 0.00 0.22mm 之間。 圖1 計(jì)算河段河勢(shì)及網(wǎng)格布置圖 Fig.1 The calculated river reach and the grids 計(jì)算區(qū)域的選取及網(wǎng)格劃分：選取塔市驛至沙夾邊長(zhǎng)約 20km 的河段作為計(jì)算河段。采用三角形網(wǎng)格，將整個(gè)計(jì)算河段劃分為 100 個(gè)斷面， 6732 個(gè)三角形網(wǎng)格單元，共有 3500 個(gè)節(jié)點(diǎn)，見圖 1 。 地形條件：以1993年10月實(shí)測(cè)河道地形為計(jì)算起始地形，并用1996年10月實(shí)測(cè)地形

12、進(jìn)行驗(yàn)證。 計(jì)算時(shí)段及水沙條件：選用1993年10月1996年10月水沙資料，按流量共劃分為120個(gè)計(jì)算時(shí)段。進(jìn)口按姚圻腦水文站流量和分組含沙量給出，出口水位根據(jù)沙夾邊與姚圻腦水位相關(guān)關(guān)系給出。 計(jì)算結(jié)果 分析 ： 圖 2 為姚圻腦計(jì)算水位與實(shí)測(cè)水位對(duì)比圖，由圖中可見計(jì)算水位與實(shí)測(cè)水位吻合較好。 圖3為不同流量下流場(chǎng)示意圖，由圖中可見大水漫灘，小水歸槽，水流平順，流場(chǎng)合理。 圖4為1996年10月計(jì)算和實(shí)測(cè)的河道地形對(duì)比圖，由圖中可見通過計(jì)算得到的地形等高線與實(shí)測(cè)河道地形等高線的位置和范圍基本吻合。 圖2計(jì)算與實(shí)測(cè)水位對(duì)比圖 Fig.2 Comparison between calculate

13、d andmeasured water level     圖3不同流量下流場(chǎng)示意圖 Fig.3 Flow fields for different discharges 圖4計(jì)算與實(shí)測(cè)河道地形比較圖 Fig.4 Comparison between calculated and measured channel     從驗(yàn)證時(shí)段的沖淤總量來看，河段內(nèi)實(shí)際淤積總量為2820萬m3，計(jì)算淤積總量為2852萬m3，相對(duì)誤差僅為1.1%。由此可見，本模型能較好反映本河段的河床變形情況。 6結(jié)論 (1) 對(duì)于有限單元法在河道平面二維水流泥沙計(jì)算中的 應(yīng)用 作了一些成功的探討，建立了一套行之有效的解法。有限元法具有能夠很好地模擬復(fù)雜邊界和河道地形等優(yōu)點(diǎn)，但若不經(jīng)過處理，就需要大量的存儲(chǔ)單元和很長(zhǎng)的運(yùn)算時(shí)間。為克服這些缺點(diǎn)，作者采用了質(zhì)量集中的方法，并提出了壓縮存儲(chǔ)的方法，使得計(jì)算存儲(chǔ)量大為減少。另外，還采用了 “ 預(yù)報(bào) - 校正 - 迭代 ” 的時(shí)間推進(jìn)算法，并提出了 “ 非恒定 - 恒定 - 非恒定流 ” 的算法，大大減少了計(jì)算量。 (2) 針對(duì)有限元前期單元剖分