黃河下游平面二維水沙運動模擬的有限元方法

1、黃河下游平面二維水沙運動模擬的有限元方法黃河下游平面二維水沙運動模擬的有限元方法摘要：本文以水流和泥沙運動 規(guī)律 的 研究 成果為基礎(chǔ)，建立了平面二維水沙模型。用有限元 方法 導(dǎo)出了本模型的離散方程式，用質(zhì)量集中和預(yù)估校正法處理、迭代求解方程組。以黃河濟南河段 1976 年汛期的洪水演進(jìn)和河床演變?yōu)槔龔乃贿^程、流量過程、斷面平均流速和最大流速、全域流速場和河床斷面沖淤形態(tài)等方面，對模型進(jìn)行了驗證模擬 計算 。數(shù)值模擬計算結(jié)果與實測或物理模型試驗結(jié)果符合較好。從而證明了本模型可靠性。關(guān)鍵詞：黃河下游 有限元方法 驗證1 前言對于河道上修建橋涵等跨河工程、引onclick=g(水工);水工程等建

2、筑物后，人們所關(guān)心的河勢變化、流速分布、河床局部沖淤形態(tài)和壅水等 問題 ，一維模型是無能為力的，但可以用二維模型解決。在進(jìn)行黃河水沙運動的模擬過程中，一些模型往往需要對方程組中的許多參數(shù)進(jìn)行經(jīng)驗處理，而且在處理復(fù)雜的河道邊界以及方程的離散和求解等方面還有許多問題值得研究。本模型以黃河onclick=g(水利);水利科研究院在泥沙運動規(guī)律方面的研究成果為基礎(chǔ)，建立模型框架。并選用適用于渾水的群體沉速公式計算泥沙沉速，引入適應(yīng)于從清水到高含沙的水流挾沙能力計算公式和動床阻力計算公式，克服數(shù)學(xué)模型參數(shù)過多而不能通用的缺陷。對模型研究區(qū)域的離散和方程的求解是數(shù)值模擬的關(guān)鍵。在離散方法方面，各種離散方法

3、都有其優(yōu)缺點，而且對每種方法又可分為多種形式，根據(jù)各種離散方法的特點和河道形態(tài)，為使對區(qū)域的離散既能很好地擬合長寬比很大，而且彎曲復(fù)雜的河道邊界，又能根據(jù)河勢、主流和水深的變化對局部區(qū)域加密細(xì)劃，本模型選用了有限元法。由于離散后形成的方程組的計算量巨大，用一般計算方法對線性方程組求解不能滿足要求，本模型 應(yīng)用 質(zhì)量集中的方法劃系數(shù)矩陣為三對角矩陣，并用預(yù)估校正法處理、迭代求解方程組。大大減少了計算量。2 基本方程和定解條件2.12.1 基本方程基本方程水流運動的基本方程為(1)(2)式中 Ui為垂線平均流速；H 為水深；Z 為水位；C 為謝才系數(shù)。 C=1/nR1/6(水力半徑 RH)，g，分

4、別為重力加速度，水密度和粘滯系數(shù)。f 為科氏力系數(shù)(f=2sin，為地球自轉(zhuǎn)角速度，為地理緯度)。為系數(shù)矩陣泥沙運動方程為1(3)(4)其中 U 為流速；i為泥沙沉速；S 和 S*分別為水流含沙量和挾沙力，f1為泥沙非飽和系數(shù)；K1為考慮紊流脈動在水平方向產(chǎn)生的擴散作用及泥沙存在產(chǎn)生的附加 影響 而引入的修正系數(shù)，簡稱為附加系數(shù)； *為平衡含沙量分布系數(shù)，詳見1。水流挾沙力計算，采用2中的公式計算水流挾沙力。河床糙率計算，應(yīng)用3中的糙率計算公式，可以描述水力泥沙因子的變化對摩阻特性的影響。2.22.2 定解條件定解條件邊界條件。對入流邊界，給出水流流速(或單寬流量)和含沙量過程；對出口邊界，給

5、出水位過程線或流速(單寬流量)過程線；對固壁邊界，法向流速為0，水流沿切線方向流速非 0。初始條件。給出在計算的初始時刻地形、流速、水位和含沙量等物理量的初始值。3 有限元離散模式的建立和方程求解3.13.1 有限元離散模式的建立有限元離散模式的建立有限元網(wǎng)格的選擇。根據(jù)需要，本次選用三角形常應(yīng)變單元類型離散 研究區(qū)域。設(shè)整個區(qū)域共劃分為 NE 個單元，I 個節(jié)點，單元節(jié)點總體編號為 i, i=1,2,I。對水沙方程(1) (4)，用 Galerkin 加權(quán)余量法逼近；對研究區(qū)域進(jìn)行剖分；對單元節(jié)點和整體節(jié)點分別編號，建立局部節(jié)點編號系統(tǒng)和整體節(jié)點編號系統(tǒng)，并確定兩個編號系統(tǒng)的關(guān)系；在離散區(qū)域

6、的基礎(chǔ)上，求出未知量在每個單元上的形函數(shù)；把形函數(shù)代入 Galerkin 積分表達(dá)式進(jìn)行單元 分析 ，建立局部有限元方程式；對所有局部有限元求和、總體合成建立總體有限元方程式，加上本質(zhì)邊界條件，即可得到本 問題 的有限元方程式。對問題進(jìn)行有限元分析式H，U，V，S 表示求解變量的變分，表示 計算 域。對以上二階導(dǎo)數(shù)項利用 Green 公式分部積分公式，并設(shè)任一單元第 i 節(jié)點的平均流速，水位和含沙量及對應(yīng)的加權(quán)函數(shù)的形函數(shù)值分別為 Ui，Vi，Hi，Si和 U*i,V*i,H*i,S*i。則有將上述表達(dá)式代入方程可得有限元方程式MijdZsi/dt=-Pij(HU)j-Pzij(HV)jMij

7、dUj/dt=-NijUj-gP1ijZsj-Mij/-(Qij+Rij)UjMijdVj/dt=-NijVj-gP2ijZsj-Mijyi/-(Qij+Rij)VjMijdSj/dt=-NijSj-s(Qij+Rij)Sj-K1*Mij(f1S-S*)/Hj i, j=1, 2, 3,I其中i,j,k=1,2,I其中對以上各式，當(dāng)在同一項中含有兩個相同的下標(biāo)時，就意味著該項表示在單元體內(nèi)全體同節(jié)點編號項的疊加；根據(jù)以上積分式先對各個單元分析，再疊加全域各單元方程式并使之滿足邊界條件即得到整體有限元方程式。對上述有限元方程在時間上對變量 Z、U、V 和 S 用向前差分格式離散可得到一個 II

8、常微分方程組。用常規(guī) 方法 對 II 階方程組求解，計算量大，本文用質(zhì)量集中方法4化方程組的系數(shù)矩陣成對角矩陣，方程組就可以容易解出。3.23.2 方程的求解方程的求解用于預(yù)估校正法迭代求解方程組第一時段的計算采用歐拉格式，第二時段采用預(yù)估校正格式，即對函數(shù)值 fn進(jìn)行預(yù)估、校正。用 fn-f*n進(jìn)行判斷，若上式成立，則 fn=f*n，否則令 fn(fn+fn-1)/2，再由上式求 f*n進(jìn)行繼續(xù)迭代，直到滿足精度要求。當(dāng)求出 Ui, Vi, Hi和 Si后代入有關(guān)方程，即求出沖淤變形的河床高程，完成方程的求解。4 模型的驗證河段基本概況 黃河下游北展滯洪區(qū)以南的北店子至后張莊河段，長約30k

9、m，是受工程控制的彎曲性河道，河道縱比降約為 1，其中北店子至濼口鐵橋河段長 20km，堤距一般 0.7km 至 1.5km。濼口鐵橋至后張莊長約 10km，堤距 1.5km2.5km，汛期平均懸移質(zhì)中值粒徑 d50=0.02mm0.027mm，河床質(zhì)中值粒徑 D50=0.07mm0.11mm，河床糙率 n=0.0130.016。模型 計算 區(qū)域的選擇和網(wǎng)格劃分。選薩口斷面為距濼口水文站上游約 20km 的北店子處，進(jìn)口斷面距其下游的曹家圈斷面、鄭家店斷面分別約 3km 和 9km；出口斷面選擇在距濼口斷面下游約 10km 的后張莊斷面附近。計算區(qū)域內(nèi)共有曹家圈、鄭家店、洛口和后張莊四個大斷面

10、，兩岸邊界選擇在險工、護(hù)灘控導(dǎo)工程和大堤等較穩(wěn)固的工程所連接的邊線上。整個計算區(qū)域劃分成 128 個小斷面，3728 個三角網(wǎng)格單元體，共有 2000 個節(jié)點，三角形沿河寬方向最小長度為 30m。一般沿河寬方向網(wǎng)格邊長為主槽 40m50m 左右，灘地 100m 左右。網(wǎng)格圖如圖 1 所示。驗證時段及水沙條件 選用 1976 年汛期 8 月 2 日至 10 月 15 日及對應(yīng)的水沙過程，共 75 天。初始條件 假設(shè)初始地形、流場、水位和含沙量為一定的常數(shù)。邊界條件 給出入流斷面各結(jié)點的流速(單寬流量)和含沙量過程和出口斷斷面的水位過程線；對兩岸邊界，根據(jù)需要分別設(shè)定滑動和不滑動邊界條件。圖 1

11、計算河段網(wǎng)格圖grids of calculated river reach驗證結(jié)果及其 分析 。圖 2 給出了鄭家店斷面和濼口斷面的水位隨時間的變化過程。從圖中可以看出，數(shù)模計算結(jié)果與實測值符合較好。圖 3 給出了濼口斷面流量過程與實測流量過程的比較圖，可以看出，計算值與實測值基本一致。圖 2 水位比較圖Water stages comparisons圖 3 流量比較圖Comparison of discharge表 1 計算流速與實測流速比較表Compasion of calculated and field velocities日期8.118.239.49.69.149.2810.2流量

12、(m3/s)1550351066107410524038502120平均流速實測2.062.622.642.832.662.412.08(m/s)計算1.862.502.502.692.532.332.00最大流速實測3.043.613.843.863.653.262.86(m/s)計算2.763.403.563.723.603.192.70圖 4 不同流量流場圖Flow fields with different dischsrge表 1 給出了一些流量的濼口斷面最大流速和最小流速的計算與實測值的比較，從圖 4 給出了流量為 2600m3/s 和 7800m3/s 下的流場的計算結(jié)果; 圖

13、4 給出了汛前小水、洪峰期大水全域流速分布，從這些圖中可以看出，不同流量下河勢、主流線的變化、水流漫灘、流速在主槽和灘地的不同分布情況。圖 5 橫斷面沖淤變形比較圖Cross-sectional erosion and deposition comparisons圖 5 給出了鄭家店、濼口重要斷面的 1976 年汛前 6 月 10 日實測地形、汛后 10 月 15 日數(shù)學(xué)模型計算地形與實測地形及物理模型試驗的對照結(jié)果。從圖中可以看出，數(shù)模計算結(jié)果在其沖淤部位和沖淤面積等方面基本上與實測數(shù)據(jù)或物理模型試驗數(shù)據(jù)相符合。5 結(jié)語1. 模型采用了黃河泥沙 研究 的新成果， 理論 基礎(chǔ)可靠。2. 采用有限元法對區(qū)域進(jìn)行離散，容易處理和更好地擬合不規(guī)則河道邊界；網(wǎng)格劃分靈活，可以對局部區(qū)域任意加密；有限元程序模塊可移植性強，很容易進(jìn)行不同網(wǎng)格形狀單元離散模式之間轉(zhuǎn)換。3. 驗證結(jié)果表明，在水位、流量、流速和河床沖淤等方面，數(shù)模計算結(jié)果與實測數(shù)據(jù)或物理模型試驗結(jié)果符合較好，從而證明了模型的可靠性。參參 考考 文文 獻(xiàn)獻(xiàn)1 張紅武，呂