水利工程論文-黃河下游平面二維水沙運動模擬的有限元方法摘要：本文以水流和泥沙運動規(guī)律的研究成果為基礎(chǔ)，建立了平面二維水沙模型。用有限元方法導(dǎo)出了本模型的離散方程式，用質(zhì)量集中和預(yù)估校正法處理、迭代求解方程組。以黃河濟南河段1976年汛期的洪水演進和河床演變?yōu)槔龔乃贿^程、流量過程、斷面平均流速和最大流速、全域流速場和河床斷面沖淤形態(tài)等方面，對模型進行了驗證模擬計算。數(shù)值模擬計算結(jié)果與實測或物理模型試驗結(jié)果符合較好。從而證明了本模型可靠性。關(guān)鍵詞：黃河下游有限元方法驗證1前言對于河道上修建橋涵等跨河工程、引水工程等建筑物后，人們所關(guān)心的河勢變化、流速分布、河床局部沖淤形態(tài)和壅水等問題，一維模型是無能為力的，但可以用二維模型解決。在進行黃河水沙運動的模擬過程中，一些模型往往需要對方程組中的許多參數(shù)進行經(jīng)驗處理，而且在處理復(fù)雜的河道邊界以及方程的離散和求解等方面還有許多問題值得研究。本模型以黃河水利科研究院在泥沙運動規(guī)律方面的研究成果為基礎(chǔ)，建立模型框架。并選用適用于渾水的群體沉速公式計算泥沙沉速，引入適應(yīng)于從清水到高含沙的水流挾沙能力計算公式和動床阻力計算公式，克服數(shù)學(xué)模型參數(shù)過多而不能通用的缺陷。對模型研究區(qū)域的離散和方程的求解是數(shù)值模擬的關(guān)鍵。在離散方法方面，各種離散方法都有其優(yōu)缺點，而且對每種方法又可分為多種形式，根據(jù)各種離散方法的特點和河道形態(tài)，為使對區(qū)域的離散既能很好地擬合長寬比很大，而且彎曲復(fù)雜的河道邊界，又能根據(jù)河勢、主流和水深的變化對局部區(qū)域加密細劃，本模型選用了有限元法。由于離散后形成的方程組的計算量巨大，用一般計算方法對線性方程組求解不能滿足要求，本模型應(yīng)用質(zhì)量集中的方法劃系數(shù)矩陣為三對角矩陣，并用預(yù)估校正法處理、迭代求解方程組。大大減少了計算量。2基本方程和定解條件2.1基本方程水流運動的基本方程為(1)(2)式中Ui為垂線平均流速；H為水深；Z為水位；C為謝才系數(shù)。C=1/nR1/6(水力半徑RH)，g，分別為重力加速度，水密度和粘滯系數(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.2定解條件邊界條件。對入流邊界，給出水流流速(或單寬流量)和含沙量過程；對出口邊界，給出水位過程線或流速(單寬流量)過程線；對固壁邊界，法向流速為0，水流沿切線方向流速非0。初始條件。給出在計算的初始時刻地形、流速、水位和含沙量等物理量的初始值。3有限元離散模式的建立和方程求解3.1有限元離散模式的建立有限元網(wǎng)格的選擇。根據(jù)需要，本次選用三角形常應(yīng)變單元類型離散研究區(qū)域。設(shè)整個區(qū)域共劃分為NE個單元，I個節(jié)點，單元節(jié)點總體編號為i,i=1,2,I。對水沙方程(1)(4)，用Galerkin加權(quán)余量法逼近；對研究區(qū)域進行剖分；對單元節(jié)點和整體節(jié)點分別編號，建立局部節(jié)點編號系統(tǒng)和整體節(jié)點編號系統(tǒng)，并確定兩個編號系統(tǒng)的關(guān)系；在離散區(qū)域的基礎(chǔ)上，求出未知量在每個單元上的形函數(shù)；把形函數(shù)代入Galerkin積分表達式進行單元分析，建立局部有限元方程式；對所有局部有限元求和、總體合成建立總體有限元方程式，加上本質(zhì)邊界條件，即可得到本問題的有限元方程式。對問題進行有限元分析式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。則有將上述表達式代入方程可得有限元方程式MijdZsi/dt=-Pij(HU)j-Pzij(HV)jMijdUj/dt=-NijUj-gP1ijZsj-Mij/-(Qij+Rij)UjMijdVj/dt=-NijVj-gP2ijZsj-Mijyi/-(Qij+Rij)VjMijdSj/dt=-NijSj-s(Qij+Rij)Sj-K1*Mij(f1S-S*)/Hji,j=1,2,3,I其中i,j,k=1,2,I其中對以上各式，當(dāng)在同一項中含有兩個相同的下標(biāo)時，就意味著該項表示在單元體內(nèi)全體同節(jié)點編號項的疊加；根據(jù)以上積分式先對各個單元分析，再疊加全域各單元方程式并使之滿足邊界條件即得到整體有限元方程式。對上述有限元方程在時間上對變量Z、U、V和S用向前差分格式離散可得到一個II常微分方程組。用常規(guī)方法對II階方程組求解，計算量大，本文用質(zhì)量集中方法4化方程組的系數(shù)矩陣成對角矩陣，方程組就可以容易解出。3.2方程的求解用于預(yù)估校正法迭代求解方程組第一時段的計算采用歐拉格式，第二時段采用預(yù)估校正格式，即對函數(shù)值fn進行預(yù)估、校正。用fn-f*n進行判斷，若上式成立，則fn=f*n，否則令fn(fn+fn-1)/2，再由上式求f*n進行繼續(xù)迭代，直到滿足精度要求。當(dāng)求出Ui,Vi,Hi和Si后代入有關(guān)方程，即求出沖淤變形的河床高程，完成方程的求解。4模型的驗證河段基本概況黃河下游北展滯洪區(qū)以南的北店子至后張莊河段，長約30km，是受工程控制的彎曲性河道，河道縱比降約為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的北店子處，進口斷面距其下游的曹家圈斷面、鄭家店斷面分別約3km和9km；出口斷面選擇在距濼口斷面下游約10km的后張莊斷面附近。計算區(qū)域內(nèi)共有曹家圈、鄭家店、洛口和后張莊四個大斷面，兩岸邊界選擇在險工、護灘控導(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é)點的流速(單寬流量)和含沙量過程和出口圖1計算河段網(wǎng)格圖gridsofcalculatedriverreach斷斷面的水位過程線；對兩岸邊界，根據(jù)需要分別設(shè)定滑動和不滑動邊界條件。驗證結(jié)果及其分析。圖2給出了鄭家店斷面和濼口斷面的水位隨時間的變化過程。從圖中可以看出，數(shù)模計算結(jié)果與實測值符合較好。圖3給出了濼口斷面流量過程與實測流量過程的比較圖，可以看出，計算值與實測值基本一致。圖2水位比較圖Waterstagescomparisons圖3流量比較圖Comparisonofdischarge表1計算流速與實測流速比較表Compasionofcalculatedandfieldvelocities日期8.11