水利工程論文-剖面二維非恒定懸移質(zhì)泥沙擴(kuò)散方程的數(shù)值方法.doc

水利工程論文-剖面二維非恒定懸移質(zhì)泥沙擴(kuò)散方程的數(shù)值方法摘要：通過(guò)討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C格式)并通過(guò)一個(gè)具體的數(shù)值例子說(shuō)明了計(jì)算的方法步驟。關(guān)鍵詞：擴(kuò)散方程差分格式精度穩(wěn)定性1引言數(shù)學(xué)模擬方法正在成為研究河流泥沙問(wèn)題的重要手段。目前，一維數(shù)學(xué)模型發(fā)展較成熟，已廣泛應(yīng)用于模擬長(zhǎng)河段的長(zhǎng)期變形，但它只能給出河段平均沖淤深度的沿程變化，如需了解短河段的河床變形細(xì)節(jié)，則要采用二維以至三維數(shù)學(xué)模型。不論是一維數(shù)學(xué)模型還是平面二維維數(shù)學(xué)模型，都不能反映含沙量沿垂線的分布狀況，并忽略了含沙量沿垂線分布對(duì)垂線平均含沙量變化過(guò)程的影響。要解決這類問(wèn)題，必須建立剖面二維數(shù)學(xué)模型。這種模型主要通過(guò)解剖面二維泥沙擴(kuò)散方程來(lái)研究懸移質(zhì)泥沙沿水深的分布及含沙量的變化過(guò)程，對(duì)水電站進(jìn)口和其它引水工程的引水口高程的確定都能提供較好的數(shù)值模擬。泥沙擴(kuò)散方程實(shí)際上是一個(gè)變系數(shù)的二階線性偏微分方程，這樣的方程在各種復(fù)雜邊界條件下求解是極為困難的。求擴(kuò)散方程的解析解在數(shù)學(xué)上存在著難以克服的困難，往往只能通過(guò)對(duì)方程的簡(jiǎn)化，才能得到一些簡(jiǎn)單邊界條件下的解析解，在這方面，A.A.Kalinske、野滿隆治、W.E.Dobbins、俞維強(qiáng)、張啟舜、韋直林等都做了有益的嘗試；求擴(kuò)散方程的數(shù)值解曾經(jīng)因?yàn)槿狈Ω咝实挠?jì)算工具而難以實(shí)現(xiàn)，直到60年代后，隨著計(jì)算機(jī)的廣泛應(yīng)用，在各種復(fù)雜邊界條件下求擴(kuò)散方程的數(shù)值解不但成為可能，而且得到迅速的發(fā)展，在這方面，曹志先、崔俠等做了大量工作，取得了很多成果。數(shù)值方法相對(duì)于解析方法在求解偏微分方程上有著明顯的優(yōu)勢(shì)，即簡(jiǎn)單靈活、計(jì)算方便快捷，但要尋找一種精度高、穩(wěn)定性好、計(jì)算方便的差分格式也并非易事。本文擬在前人研究的基礎(chǔ)上著重討論剖面二維泥沙擴(kuò)散方程的數(shù)值解問(wèn)題，希望能提供一種精度高、穩(wěn)定性好、計(jì)算方便的數(shù)值解。2基本方程剖面二維泥沙擴(kuò)散方程的形式為(1)式中x,y為水流方向和鉛直方向的維軸；u,v分別為沿水流方向和鉛直方向的時(shí)均流速；sx,sy分別為水流方向和鉛直方向的泥沙擴(kuò)散系數(shù)；,S分別為泥沙靜水沉速和含沙量。對(duì)于式(1)的求解，研究者一般會(huì)對(duì)它進(jìn)行不同程度的簡(jiǎn)化，為此引入以下假定中的一種或幾種：A.非恒定流可以概化為梯級(jí)式恒定流，即；B.在一個(gè)時(shí)段內(nèi)，認(rèn)為泥沙運(yùn)動(dòng)可以概化為處于恒定狀態(tài)，即；C.在二維流動(dòng)中，縱向擴(kuò)散系數(shù)與方程其他項(xiàng)相比，可以忽略不計(jì)，即認(rèn)為方程右端第一項(xiàng)可以忽略；D.認(rèn)為懸移質(zhì)泥沙粒徑均一，即=const；E.認(rèn)為水流為二維均勻流，即v=0。為簡(jiǎn)單起見(jiàn)，我們討論的范圍限于水流條件為二維非恒定均勻流，懸移質(zhì)泥沙粒徑均勻，為此引入假定C、D、E。這時(shí)，泥沙擴(kuò)散方程為(2)目前，對(duì)s的變化規(guī)律研究得不很充分，一般假定s=m(3)其中m為動(dòng)量傳遞系數(shù)，為修正值。由勃蘭特爾摻長(zhǎng)理論可得m=u*y/(1-y/h)(4)式中為卡門常數(shù)，u*為摩阻流速。對(duì)于u，我們?nèi)】?勃蘭特爾對(duì)數(shù)流速分布公式(umax-u)/u*=(1/)ln(h/y)(5)令W=+(u*/h)(1-2y/h),(2)可變形為(6)3差分方程3.1網(wǎng)格的剖分t=，X方向的x=h1，Y方向的空間步長(zhǎng)為y=h2,這樣形成如下網(wǎng)格Dh=(xj,yl,tn)xj=jh1,yl=lh2,tn=n其中j=0,.,N；l=0,.,M；n0；3.2構(gòu)造差分格式通過(guò)對(duì)流方程和擴(kuò)散方程的差分格式的構(gòu)造，我們可以得到對(duì)流擴(kuò)散方程的差分格式。由于隱式格式穩(wěn)定性好，考慮Crank-Nicholson型隱式格式。為此，引入差分算子記號(hào)2ySnj,l=Snj,l+1-2Snj,l+Snj,l-1LxSnj,l=Snj+1,l-Snj-1,lLySnj,l=Snj,l+1-Snj,l-1為了看得更清楚，暫且取h1=h2=h.對(duì)式(6)離散，則C-N格式為Sn+1j,l-Snj,l/+u/4hLx(Snj,l+Sn+1j,l)=s/2h22y(Sn+1j,l+Snj,l)+W/4hLy(Sn+1j,l+Snj,l)(7)C-N格式的精度是二階的，絕對(duì)穩(wěn)定。但對(duì)于二維問(wèn)題，由(7)導(dǎo)出的方程組，其系數(shù)矩陣不是三對(duì)角矩陣，不能用追趕法求解。因此，考慮構(gòu)造交替方向的隱式格式(命名為Z-C格式)(Sn+1/2j,l-Snj,l)/(/2)+u/2hLxSnj,l=s/h22ySn+1/2j,l+W/2hLySn+1/2j,l(8)(Sn+1j,l-Sn+1/2j,l)/2+u/2hLxSn+1j,l=s/h22ySn+1/2j,l+W/2hLySn+1/2j,l(9)可以看出，計(jì)算Sn+1j,l是由兩步組成的，每一步僅是一個(gè)方向的隱式，故用兩次追趕法即可。3.3精度分析現(xiàn)在，我們考慮Z-C格式的精度。先設(shè)法消去過(guò)渡值Sn+1/2j,l，為此，將(8)和(9)兩式相加，可得Sn+1j,l-Snj,l+u/2hLx(Sn+1j,l+Snj,l)=s/2h22ySn+1/2j,l+W/4hLxSn+1/2j,l(10)將(8)和(9)兩式相減，可得Sn+1/2j,l=(Sn+1j,l+Snj,l)/2+s/4h22y(Snj,l-Sn+1j,l)+W/8hLy(Snj,l-Sn+1j,l)(11)把式(11)代入(10)，變形整理，可得(Sn+1j,l-Snj,l)(1-us/8h3Lx2y-2uW/16h2LxLy)=(Snj,l+Sn+1j,l)(us/2h22y+W/4hLy-u/4hLx)(12)設(shè)S(x,y,t)是(12)的精確解，并假定S(x,y,t)關(guān)于t三次連續(xù)可微，關(guān)于x,y四次連續(xù)可微，那么利用Taylor級(jí)數(shù)展開可得S(xj,yl,tn+1)-S(xj,yl,tn)/(1-2s/8h3Lx2y-2uW/16h2LxLy)-S(xj,yl,tn+1)-S(xj,yl,tn)/(s/2h2)2y+W/4hLy-u/4hLx)=O(2+h2)(13)由此可見(jiàn)，Z-C格式具有二階精度。3.4穩(wěn)定性分析現(xiàn)在，我們來(lái)討論Z-C格式的穩(wěn)定性。為此，把式(12)變形整理得(1-s/2h22y-W/4hLy)(1+u/4hLx)Sn+1j,l=(1+s/2h2y2+W/4hLy)(1-u/4hLx)Snj,l(14)由式(14)可得出過(guò)渡因子為G(,k)=(1-2s/h2sin2k2h/2+iW/2hsink2h)(1-iu/2hsink1h)/(1+2s/h2sin2k2h/2-iW/2hsink2h)(1+iu/2hsink1h)(15)令a=2s/h2sin2k2h/2,b=W/2hsink2h,c=u/2hsink1h,則|G|2=|1-a2-b2+2bi/(1+a)2+b2|2|1-c2-2ci/1+c2|2=1(1-a2-b2)2+4b2/(1+a)2+b22=1-4a2+4a+4b2+4a3/(1+a)2+b22(16)顯然，對(duì)于任意的,h,|G(,k)|21,Z-C格式是絕對(duì)穩(wěn)定的。4數(shù)值計(jì)算4.1邊界條件我們考慮初邊值問(wèn)題。(1)初始條件用Rouse公式給出含沙量沿垂線分布S(x,y,0)=Sa5(a/h-a)z(h-y/y)z(17)式中z=/ku*為懸浮指標(biāo)，Sa為近底含沙量，h為水深，一般取a=0.010.05h。(2)水面條件(y=h)(18)(3)底部邊界條件(y=a)(19)式中Sa*為近底挾沙力，即輸沙平衡時(shí)的近底售沙量Sa.(4)近底含沙量計(jì)算近底含沙量在求解泥沙擴(kuò)散方程時(shí)具有邊界條件性質(zhì)，它選取的正確與否，意味著所給邊界條件是否正確。實(shí)際工程中一般缺乏實(shí)測(cè)資料，近底含沙量不易測(cè)定。這里，我們利用水流挾沙力和含沙量沿垂線分布公式來(lái)反求近底含沙量。已知斷面平均挾沙力為S*=k(u3/gh)m(20)假定Sa*=S*(21)輸沙平衡時(shí)，含沙量沿垂線分布用Rouse公式(17)表示，用(17)表達(dá)挾沙力的垂線分布，然后沿垂線積分得斷面平均挾沙力為(22)將(22)與(21)比較，可得(23)4.2計(jì)算步驟為方便計(jì)算，將式(8)和(9)式變形整理，并對(duì)X，Y方向取不同的空間步長(zhǎng)(s/2h22-W/4h2)Sn+1/2j,l+1-(s/h22+1)Sn+1/2j,l+(s/2h22-W/4h2)Sn+1/2j,l+1=u/4h1LxSnj,l-Snj,l(24)-u/4h1Snj-1,l+Snj,l+u/4h1Snj+1,l=s/2h222ySn+1/2j,l+W/4h2LySn+1/2j,l+Sn+1/2j,l(25)在一個(gè)時(shí)間層(第n層)內(nèi)，計(jì)算分兩步進(jìn)行：第一步，對(duì)式(24)用追趕法求第n+1/2層的過(guò)渡值。令C1=-u/4h1,C2=1,C3=u/4h1,E1=s/2h22,E2=W/4h2D1=s/2h22-W/4h2,D2=-s/h22-1,D3=s/2h22+W/4h2l=1時(shí)，D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1l=2時(shí)，D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1l=M時(shí),D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3LxSnj,M-Snj,M令Hl=-Snj,l+C3LxSnj,l(1lM)H1=-Snj,1+C3LxSnj,1-D1Sn+1/2j,0HM=-Snj,M+C3LxSnj,M-D3Sn+1/2j,M+1其中Sn+1/2j,0和Sn+1/2j,M+1由邊界條件給出，則用矩陣形式表示為（26）第二步，再對(duì)(25)式用追趕法求第n+1層的值令Fj=E12ySn+1/2j,l+E2LySn+1/2j,l+Sn+1/2j,l(1jN)F1=E12ySn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C1Sn+10,lFN=E12ySn+1/2N,l+E2LySn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l其中Sn+10,l和Sn+1N+1，l由邊界條件給出，則同理可得矩陣方程(27)這樣，按此步驟一層層地計(jì)算。4.3數(shù)值模擬合理性分析受所掌握的實(shí)測(cè)資料的限制，目前尚無(wú)法對(duì)本文提出的算法與含沙量沿垂線分布的實(shí)測(cè)值進(jìn)行對(duì)比。我們用庫(kù)里阿雷克沉沙池的實(shí)測(cè)資料作了垂線平均值沿程變化的比較。該沉沙池的主要數(shù)據(jù)為：池深h=1.53m；平均流速u=0.1