水利工程論文-黃河下游平面二維水沙運動模擬的有限元方法摘要：本文以水流和泥沙運動規(guī)律的研究成果為基礎，建立了平面二維水沙模型。用有限元方法導出了本模型的離散方程式，用質量集中和預估校正法處理、迭代求解方程組。以黃河濟南河段1976年汛期的洪水演進和河床演變?yōu)槔龔乃贿^程、流量過程、斷面平均流速和最大流速、全域流速場和河床斷面沖淤形態(tài)等方面，對模型進行了驗證模擬計算。數值模擬計算結果與實測或物理模型試驗結果符合較好。從而證明了本模型可靠性。關鍵詞：黃河下游有限元方法驗證1前言對于河道上修建橋涵等跨河工程、引水工程等建筑物后，人們所關心的河勢變化、流速分布、河床局部沖淤形態(tài)和壅水等問題，一維模型是無能為力的，但可以用二維模型解決。在進行黃河水沙運動的模擬過程中，一些模型往往需要對方程組中的許多參數進行經驗處理，而且在處理復雜的河道邊界以及方程的離散和求解等方面還有許多問題值得研究。本模型以黃河水利科研究院在泥沙運動規(guī)律方面的研究成果為基礎，建立模型框架。并選用適用于渾水的群體沉速公式計算泥沙沉速，引入適應于從清水到高含沙的水流挾沙能力計算公式和動床阻力計算公式，克服數學模型參數過多而不能通用的缺陷。對模型研究區(qū)域的離散和方程的求解是數值模擬的關鍵。在離散方法方面，各種離散方法都有其優(yōu)缺點，而且對每種方法又可分為多種形式，根據各種離散方法的特點和河道形態(tài)，為使對區(qū)域的離散既能很好地擬合長寬比很大，而且彎曲復雜的河道邊界，又能根據河勢、主流和水深的變化對局部區(qū)域加密細劃，本模型選用了有限元法。由于離散后形成的方程組的計算量巨大，用一般計算方法對線性方程組求解不能滿足要求，本模型應用質量集中的方法劃系數矩陣為三對角矩陣，并用預估校正法處理、迭代求解方程組。大大減少了計算量。2基本方程和定解條件2.1基本方程水流運動的基本方程為(1)(2)式中Ui為垂線平均流速；H為水深；Z為水位；C為謝才系數。C=1/nR1/6(水力半徑RH)，g，分別為重力加速度，水密度和粘滯系數。f為科氏力系數(f=2sin，為地球自轉角速度，為地理緯度)。為系數矩陣泥沙運動方程為1(3)(4)其中U為流速；i為泥沙沉速；S和S*分別為水流含沙量和挾沙力，f1為泥沙非飽和系數；K1為考慮紊流脈動在水平方向產生的擴散作用及泥沙存在產生的附加影響而引入的修正系數，簡稱為附加系數；*為平衡含沙量分布系數，詳見1。水流挾沙力計算，采用2中的公式計算水流挾沙力。河床糙率計算，應用3中的糙率計算公式，可以描述水力泥沙因子的變化對摩阻特性的影響。2.2定解條件邊界條件。對入流邊界，給出水流流速(或單寬流量)和含沙量過程；對出口邊界，給出水位過程線或流速(單寬流量)過程線；對固壁邊界，法向流速為0，水流沿切線方向流速非0。初始條件。給出在計算的初始時刻地形、流速、水位和含沙量等物理量的初始值。3有限元離散模式的建立和方程求解3.1有限元離散模式的建立有限元網格的選擇。根據需要，本次選用三角形常應變單元類型離散研究區(qū)域。設整個區(qū)域共劃分為NE個單元，I個節(jié)點，單元節(jié)點總體編號為i,i=1,2,I。對水沙方程(1)(4)，用Galerkin加權余量法逼近；對研究區(qū)域進行剖分；對單元節(jié)點和整體節(jié)點分別編號，建立局部節(jié)點編號系統和整體節(jié)點編號系統，并確定兩個編號系統的關系；在離散區(qū)域的基礎上，求出未知量在每個單元上的形函數；把形函數代入Galerkin積分表達式進行單元分析，建立局部有限元方程式；對所有局部有限元求和、總體合成建立總體有限元方程式，加上本質邊界條件，即可得到本問題的有限元方程式。對問題進行有限元分析式H，U，V，S表示求解變量的變分，表示計算域。對以上二階導數項利用Green公式分部積分公式，并設任一單元第i節(jié)點的平均流速，水位和含沙量及對應的加權函數的形函數值分別為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其中對以上各式，當在同一項中含有兩個相同的下標時，就意味著該項表示在單元體內全體同節(jié)點編號項的疊加；根據以上積分式先對各個單元分析，再疊加全域各單元方程式并使之滿足邊界條件即得到整體有限元方程式。對上述有限元方程在時間上對變量Z、U、V和S用向前差分格式離散可得到一個II常微分方程組。用常規(guī)方法對II階方程組求解，計算量大，本文用質量集中方法4化方程組的系數矩陣成對角矩陣，方程組就可以容易解出。3.2方程的求解用于預估校正法迭代求解方程組第一時段的計算采用歐拉格式，第二時段采用預估校正格式，即對函數值fn進行預估、校正。用fn-f*n進行判斷，若上式成立，則fn=f*n，否則令fn(fn+fn-1)/2，再由上式求f*n進行繼續(xù)迭代，直到滿足精度要求。當求出Ui,Vi,Hi和Si后代入有關方程，即求出沖淤變形的河床高程，完成方程的求解。4模型的驗證河段基本概況黃河下游北展滯洪區(qū)以南的北店子至后張莊河段，長約30km，是受工程控制的彎曲性河道，河道縱比降約為1，其中北店子至濼口鐵橋河段長20km，堤距一般0.7km至1.5km。濼口鐵橋至后張莊長約10km，堤距1.5km2.5km，汛期平均懸移質中值粒徑d50=0.02mm0.027mm，河床質中值粒徑D50=0.07mm0.11mm，河床糙率n=0.0130.016。模型計算區(qū)域的選擇和網格劃分。選薩口斷面為距濼口水文站上游約20km的北店子處，進口斷面距其下游的曹家圈斷面、鄭家店斷面分別約3km和9km；出口斷面選擇在距濼口斷面下游約10km的后張莊斷面附近。計算區(qū)域內共有曹家圈、鄭家店、洛口和后張莊四個大斷面，兩岸邊界選擇在險工、護灘控導工程和大堤等較穩(wěn)固的工程所連接的邊線上。整個計算區(qū)域劃分成128個小斷面，3728個三角網格單元體，共有2000個節(jié)點，三角形沿河寬方向最小長度為30m。一般沿河寬方向網格邊長為主槽40m50m左右，灘地100m左右。網格圖如圖1所示。驗證時段及水沙條件選用1976年汛期8月2日至10月15日及對應的水沙過程，共75天。初始條件假設初始地形、流場、水位和含沙量為一定的常數。邊界條件給出入流斷面各結點的流速(單寬流量)和含沙量過程和出口圖1計算河段網格圖gridsofcalculatedriverreach斷斷面的水位過程線；對兩岸邊界，根據需要分別設定滑動和不滑動邊界條件。驗證結果及其分析。圖2給出了鄭家店斷面和濼口斷面的水位隨時間的變化過程。從圖中可以看出，數模計算結果與實測值符合較好。圖3給出了濼口斷面流量過程與實測流量過程的比較圖，可以看出，計算值與實測值基本一致。圖2水位比較圖Waterstagescomparisons圖3流量比較圖Comparisonofdischarge表1計算流速與實測流速比較表Compasionofcalculatedandfieldvelocities日期8.11