水利工程論文-非恒定流泥沙數(shù)學(xué)模型原理及其應(yīng)用摘要：本文介紹了一維非恒定泥沙數(shù)學(xué)模型的原理和計(jì)算方法，并且對(duì)于泥沙數(shù)學(xué)模型中的一些關(guān)鍵技術(shù)等問(wèn)題進(jìn)行了較詳細(xì)的討論。利用已建立的一維非恒定泥沙數(shù)學(xué)模型，對(duì)黃河小北干流1981年洪水資料進(jìn)行了模擬計(jì)算，計(jì)算結(jié)果表明：黃河干流、渭河、北洛河計(jì)算水位、流量過(guò)程線與實(shí)測(cè)值符合良好，可用于黃河的洪水演進(jìn)計(jì)算。關(guān)鍵詞：非恒定流泥沙數(shù)學(xué)模型河床變形1泥沙數(shù)學(xué)模型基本方程明渠或天然河流常被考慮作為一維流動(dòng)，根據(jù)洪水波運(yùn)動(dòng)的圣維南方程、泥沙連續(xù)方程和泥沙擴(kuò)散方程，可以簡(jiǎn)化推導(dǎo)出一維非恒定流泥沙數(shù)學(xué)模型的基本方程，其形式為水流連續(xù)方程(1)水流動(dòng)量方程(2)泥沙連續(xù)方程(3)不平衡輸沙方程(4)挾沙力方程S*=f(Q,A,B,S)(5)其中Q為流量，A為斷面面積，B為斷面寬度，Z為水位，K為流量模數(shù)，S為斷面平均含沙量，S*為斷面平均挾沙力，g為重力加速度，為泥沙非平衡恢復(fù)飽和系數(shù)，為泥沙顆粒沉速，Ad為斷面沖淤面積，s為泥沙干容重。2計(jì)算方法簡(jiǎn)介一維非恒定泥沙數(shù)學(xué)模型的計(jì)算采用非耦合方法，首先求解水流連續(xù)方程和動(dòng)量方程，然后求解水流挾沙力、泥沙不平衡輸沙方程和泥沙連續(xù)方程，具體求解過(guò)程如下。2.1水流方程的求解首先利用Preissmann隱式差分格式將水流連續(xù)方程和動(dòng)量方程轉(zhuǎn)化為差分方程，然后對(duì)差分方程進(jìn)行線性化，在線性化過(guò)程中，略去增量的乘積項(xiàng)，最后得到以下線性方程組A1jQj+B1jZj+C1jQj+1+D1jZj+1=E1j(6)A2jQj+B2jZj+C2jQj+1+D2jZj+1=E2j(7)其中Aij、Bij、Cij、Dij、Eij（i=1,2）為第j單元河段差分方程的系數(shù)（j=1,2,.,N-1,其中N為斷面?zhèn)€數(shù)）。給定邊界條件Q1=Q1n+1-Qn1=Q1(tn+1)-Q1n(8)ZN=ZNn+1-ZNn=ZN(tn+1)-ZNn(9)方程(6)、(7)式及邊界條件，共有2N個(gè)未知數(shù)，2N個(gè)方程，可以求解。由于差分方程中的系數(shù)包含有未知數(shù)，方程求解不能直接求出未知變量，因此方程求解時(shí)必須進(jìn)行迭代處理。下面給出用追趕法求解的步驟，追趕方程為Qj=FjZj+Gj(10)Zj=HjQj+1+IjZj+1+Jj(11)其中Hj、Ij、Jj、，F(xiàn)j、Gj為追趕系數(shù)。2.2水流挾沙力、動(dòng)床阻力、河寬變化、床沙級(jí)配調(diào)整采用文獻(xiàn)2的方法2.3不平衡輸沙方程求解利用迎風(fēng)格式，將(3)式離散為差分方程，整理后得(12)當(dāng)Q0時(shí)，利用上邊界條件，自上而下計(jì)算各斷面含沙量；當(dāng)Q0時(shí)，利用下邊界條件由下至上計(jì)算各斷面含沙量。2.4河床變形及淤積量計(jì)算由式(3)與式(4)相減可以得到河床變形方程為(13)將上式寫(xiě)成差分方程，各斷面淤積面積為(14)第j河段的淤積量為Wj=(Aj+Aj+1)xj/2(15)3水流內(nèi)邊界的處理水流內(nèi)邊界是指河道的幾何形狀的不連續(xù)或水力特性的不連續(xù)點(diǎn)。例如，河流的匯合點(diǎn)、河流分流、局部河段內(nèi)生產(chǎn)堤決口等等。在這些內(nèi)部邊界處，圣維南方程組和單一河道泥沙不平衡輸沙方程等都不再適用，必須根據(jù)其水力特性作特殊處理。內(nèi)邊界條件通常包含兩個(gè)相容條件：即流量的連續(xù)條件和能量守恒條件（或動(dòng)量守恒條件）。本模型主要考慮了以下幾個(gè)類(lèi)型的內(nèi)邊界處理。3.1水沙的匯入或匯出如圖1所示，假設(shè)匯入或匯出點(diǎn)上下斷面滿足以下條件3.2支流從干流分流如圖2所示，干流和支流上斷面之間應(yīng)滿足連續(xù)方程和能量方程圖1水沙的匯入或匯出示意圖Sketchofwater/sedimentinfloworoutflow圖2支流從干流分流示意圖Outflowfromthemainstem圖3支流匯入干流示意圖Inflowintothemainstem支流從干流分流時(shí)，將干流分流斷面按干支流流量比分為兩部分，忽略時(shí)變項(xiàng)，將方程(4)式直接寫(xiě)成差分形式，求得干流和支流下游斷面的含沙量為其中為分流系數(shù)，由實(shí)測(cè)資料確定。3.3支流匯入干流當(dāng)支流匯入干流時(shí)，與支流從干流分流類(lèi)似，干流和支流上斷面之間也滿足連續(xù)方程和能量方程將匯流斷面按干支流流量比分為兩部分，忽略時(shí)變項(xiàng)，將方程(4)式直接寫(xiě)成差分形式，求得匯流斷面的含沙量為其中為分流系數(shù)，由實(shí)測(cè)資料確定。4泥沙數(shù)學(xué)模型的應(yīng)用4.1計(jì)算區(qū)域和時(shí)段計(jì)算區(qū)域?yàn)辇堥T(mén)、華縣、河津、狀頭至潼關(guān)河段，黃河干流上有渭河和汾河匯入，在渭河上有北洛河匯入。本模型同時(shí)模擬黃河干流、渭河、北洛河三條河流的洪水演進(jìn)過(guò)程，各匯入點(diǎn)作為內(nèi)邊界處理，汾河僅作為已知水、沙過(guò)程線匯NFDA1控制的已知條件，下邊界條件為潼關(guān)站出口控制水位。計(jì)算時(shí)段為1981年汛初第一場(chǎng)洪水，洪水時(shí)間為7月3號(hào)至7月14號(hào)，洪水持續(xù)時(shí)間為12天，龍門(mén)最大流量為6400m3/s、最大含沙量為298.0kg/m3；華縣最大流量為970m3/s、最大含沙量為117.0kg/m3。原始大斷面資料采用1981年汛前實(shí)測(cè)大斷面資料。4.2計(jì)算結(jié)果與實(shí)測(cè)值的比較圖4為潼關(guān)站計(jì)算流量與實(shí)測(cè)值的比較、計(jì)算含沙量與實(shí)測(cè)值的比較。由圖4可以看出，計(jì)算的潼關(guān)出口流量過(guò)程線與實(shí)測(cè)過(guò)程線比較符合，計(jì)算洪峰最大值和相應(yīng)洪峰傳播時(shí)間與實(shí)測(cè)值比較接近；計(jì)算潼關(guān)出口含沙量過(guò)程線與實(shí)測(cè)值也是比較一致的。由圖4分析得出，本模型能夠比較好地模擬出三條河流的洪水傳播過(guò)程和泥沙沖淤調(diào)整過(guò)程。圖4潼關(guān)站流量、含沙量計(jì)算值與實(shí)測(cè)值的比較ComparisonofcalculatedandmeasureddischargeandconcentrationsatTongguanstation圖5華陰站流量、水位計(jì)算值與實(shí)測(cè)值的比較ComparisonofcalculateddischargeandwaterlevelswithmeasuredatHuayinstation圖5為華陰站計(jì)算水位與實(shí)測(cè)值的比較、計(jì)算流量與實(shí)測(cè)值的比較。由圖中可以看出，計(jì)算水位、流量過(guò)程線與相應(yīng)實(shí)測(cè)過(guò)程線比較符合。同時(shí)也看出，在洪水初期，華陰站出現(xiàn)了倒灌現(xiàn)象，流量出現(xiàn)了負(fù)值，本模型也比較好地模擬出了黃