黃河下游平面二維水沙運(yùn)動(dòng)模擬的有限元方法

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

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

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

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

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

6、的基礎(chǔ)上，求出未知量在每個(gè)單元上的形函數(shù)；把形函數(shù)代入 Galerkin 積分表達(dá)式進(jìn)行單元 分析 ，建立局部有限元方程式；對(duì)所有局部有限元求和、總體合成建立總體有限元方程式，加上本質(zhì)邊界條件，即可得到本 問題 的有限元方程式。對(duì)問題進(jìn)行有限元分析式H，U，V，S 表示求解變量的變分，表示 計(jì)算 域。對(duì)以上二階導(dǎo)數(shù)項(xiàng)利用 Green 公式分部積分公式，并設(shè)任一單元第 i 節(jié)點(diǎn)的平均流速，水位和含沙量及對(duì)應(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其中對(duì)以上各式，當(dāng)在同一項(xiàng)中含有兩個(gè)相同的下標(biāo)時(shí)，就意味著該項(xiàng)表示在單元體內(nèi)全體同節(jié)點(diǎn)編號(hào)項(xiàng)的疊加；根據(jù)以上積分式先對(duì)各個(gè)單元分析，再疊加全域各單元方程式并使之滿足邊界條件即得到整體有限元方程式。對(duì)上述有限元方程在時(shí)間上對(duì)變量 Z、U、V 和 S 用向前差分格式離散可得到一個(gè) II

8、常微分方程組。用常規(guī) 方法 對(duì) II 階方程組求解，計(jì)算量大，本文用質(zhì)量集中方法4化方程組的系數(shù)矩陣成對(duì)角矩陣，方程組就可以容易解出。3.23.2 方程的求解方程的求解用于預(yù)估校正法迭代求解方程組第一時(shí)段的計(jì)算采用歐拉格式，第二時(shí)段采用預(yù)估校正格式，即對(duì)函數(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 模型的驗(yàn)證河段基本概況 黃河下游北展滯洪區(qū)以南的北店子至后張莊河段，長(zhǎng)約30k

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

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

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

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

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