剖面二維非恒定懸移質(zhì)泥沙擴(kuò)散方程的數(shù)值方法(1)

1、    剖面二維非恒定懸移質(zhì)泥沙擴(kuò)散方程的數(shù)值方法(1)    通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C格式)并通過一個具體的數(shù)值例子說明了計算的方法步驟。 關(guān)鍵詞：擴(kuò)散方程 差分格式 精度 穩(wěn)定性 1 引言數(shù)學(xué)模擬方法正在成為研究河流泥沙問題的重手段。目前，一維數(shù)學(xué)模型發(fā)展較成熟，已廣泛應(yīng)用于模擬長河段的長期變形，但它只能給出河段平均沖淤深度的沿程變化，如需了解短河段的河床變形細(xì)節(jié)，則采用二維以至三維數(shù)學(xué)模型。不論是一維數(shù)學(xué)模型還是平面二維維數(shù)學(xué)模型，都

2、不能反映含沙量沿垂線的分布狀況，并忽略了含沙量沿垂線分布對垂線平均含沙量變化過程的影響。解決這類問題，必須建立剖面二維數(shù)學(xué)模型。這種模型主通過解剖面二維泥沙擴(kuò)散方程來研究懸移質(zhì)泥沙沿水深的分布及含沙量的變化過程，對水電站進(jìn)口和其它引水工程的引水口高程的確定都能提供較好的數(shù)值模擬。泥沙擴(kuò)散方程實際上是一個變系數(shù)的二階線性偏微分方程，這樣的方程在各種復(fù)雜邊界條件下求解是極為困難的。求擴(kuò)散方程的解析解在數(shù)學(xué)上存在著難以克服的困難，往往只能通過對方程的簡化，才能得到一些簡單邊界條件下的解析解，在這方面，A.A.Kalinske、野滿隆治、W.E.Dobbins、俞維強、張啟舜、韋直林等都做了有益的嘗試

3、；求擴(kuò)散方程的數(shù)值解曾經(jīng)因為缺乏高效率的計算工具而難以實現(xiàn)，直到60年代后，隨著計算機(jī)的廣泛應(yīng)用，在各種復(fù)雜邊界條件下求擴(kuò)散方程的數(shù)值解不但成為可能，而且得到迅速的發(fā)展，在這方面，曹志先、崔俠等做了大量工作，取得了很多成果。數(shù)值方法相對于解析方法在求解偏微分方程上有著明顯的優(yōu)勢，即簡單靈活、計算方便快捷，但尋找一種精度高、穩(wěn)定性好、計算方便的差分格式也并非易事。本文擬在前人研究的基礎(chǔ)上著重討論剖面二維泥沙擴(kuò)散方程的數(shù)值解問題，希望能提供一種精度高、穩(wěn)定性好、計算方便的數(shù)值解。2 基本方程         剖面二維泥

4、沙擴(kuò)散方程的形式為            (1)式中 x,y為水流方向和鉛直方向的維軸；u,v分別為沿水流方向和鉛直方向的時均流速；sx,sy分別為水流方向和鉛直方向的泥沙擴(kuò)散系數(shù)；,S分別為泥沙靜水沉速和含沙量。對于式(1)的求解，研究者一般會對它進(jìn)行不同程度的簡化，為此引入以下假定中的一種或幾種：     A.非恒定流可以概化為梯級式恒定流，即       

5、60; B.在一個時段內(nèi)，認(rèn)為泥沙運動可以概化為處于恒定狀態(tài)，即    C.在二維流動中，縱向擴(kuò)散系數(shù)與方程其他項相比，可以忽略不計，即認(rèn)為方程右端第一項可以忽略；D.認(rèn)為懸移質(zhì)泥沙粒徑均一，即=const;E.認(rèn)為水流為二維均勻流，即v=0。為簡單起見，我們討論的范圍限于水流條件為二維非恒定均勻流，懸移質(zhì)泥沙粒徑均勻，為此引入假定C、D、E。這時，泥沙擴(kuò)散方程為         (2)    目前，對s的變化規(guī)律研究得不很充分，一般假

6、定         s=m    (3)其中m為動量傳遞系數(shù)，為修正值。由勃蘭特爾摻長理論可得     m=u*y/(1-y/h)    (4)    式中為卡門常數(shù)，u*為摩阻流速。    對于u，我們?nèi)】?勃蘭特爾對數(shù)流速分布公式     (umax-u)/u*=(1/)l

7、n(h/y)    (5)令W= (u            摘通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C         本篇論文是由3COME文檔頻道的網(wǎng)友為您在網(wǎng)絡(luò)上收集整理餅投稿至本站的，論文版權(quán)屬原作者，請不用于商業(yè)用途或者抄襲，僅供參考學(xué)習(xí)之用，否者后果自負(fù)，如果此文侵犯您的合法權(quán)益，請聯(lián)

8、系我們。            摘通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C         本篇論文是由3COME文檔頻道的網(wǎng)友為您在網(wǎng)絡(luò)上收集整理餅投稿至本站的，論文版權(quán)屬原作者，請不用于商業(yè)用途或者抄襲，僅供參考學(xué)習(xí)之用，否者后果自負(fù)，如果此文侵犯您的合法權(quán)益，請聯(lián)系我們。j,l-Snj,l)(1-us/8h3Lx2y-2uW/16h

9、2LxLy)=(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級數(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)   

10、60;(13)由此可見，Z-C格式具有二階精度。3.4 穩(wěn)定性分析現(xiàn)在，我們來討論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)可得出過渡因子為     G(,k)=(1-2s/h2sin2k2h/2 iW/2hsink2h)(1-iu/2hsink1h)/(1 2s/h2sin2k2h/2-iW/2

11、hsink2h)(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)顯然，對于任意的,h,|G(,k)|2<1,所以Z-C格式是絕對穩(wěn)定的。

12、    4 數(shù)值計算4.1 邊界條件我們考慮初邊值問題。(1)初始條件     用Rouse公式給出含沙量沿垂線分布        S(x,y,0)=Sa·5(a/h-a)z(h-y/y)z    (17)    式中z=/ku*為懸浮指標(biāo)，Sa為近底含沙量，h為水深，一般取a=0.010.05h。    &

13、#160;    (2)水面條件        (y=h)    (18)    (3)底部邊界條件         (y=a)    (19)式中Sa*為近底挾沙力，即輸沙平衡時的近底售沙量Sa.(4)近底含沙量計算近底含沙量在求解泥沙擴(kuò)散方程時具有邊界條件性質(zhì)，它選取的正確與否，意味著

14、所給邊界條件是否正確。實際工程中一般缺乏實測資料，近底含沙量不易測定。這里，我們利用水流挾沙力和含沙量沿垂線分布公式來反求近底含沙量。已知斷面平均挾沙力為     S*            摘通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C         本篇論文是由3COME文檔頻道的網(wǎng)友為您在網(wǎng)絡(luò)上

15、收集整理餅投稿至本站的，論文版權(quán)屬原作者，請不用于商業(yè)用途或者抄襲，僅供參考學(xué)習(xí)之用，否者后果自負(fù)，如果此文侵犯您的合法權(quán)益，請聯(lián)系我們。=k(u3/gh)m    (20)    假定         Sa*=S*    (21)輸沙平衡時，含沙量沿垂線分布用Rouse公式(17)表示，用(17)表達(dá)挾沙力的垂線分布 ，然后沿垂線積分得斷面平均挾沙力為   

16、0;     (22)    將(22)與(21)比較，可得             (23)4.2 計算步驟為方便計算，將式(8)和(9)式變形整理，并對X，Y方向取不同的空間步長     (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-

17、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)在一個時間層(第n層)內(nèi)，計算分兩步進(jìn)行：第一步，對式(24)用追趕法求第n 1/2層的過渡值。     令C1=-u/4h1,C2=1,C3=u/4h1,E1=s/2h22,E2=W/4h2   

18、0;D1=s/2h22-W/4h2,D2=-s/h22-1,D3=s/2h22 W/4h2    l=1時，D1Sn 1/2j,0 D2Sn 1/2j,1 D3Sn 1/2j,2=C3LxSnj,1-Snj,1    l=2時，D1Sn 1/2j,0 D2Sn 1/2j,1 D3Sn 1/2j,2=C3LxSnj,1-Snj,1        l=M時,D1Sn 1/2j,M-1 D2Sn 1/2j,M D3Sn 1/2j,M 1=C3

19、LxSnj,M-Snj,M    令Hl=-Snj,l C3LxSnj,l (1<l<M)    H1=-Snj,1 C3LxSnj,1-D1Sn 1/2j,0    HM=-Snj,M C3LxSnj,M-D3Sn 1/2j,M 1    其中Sn 1/2j,0和Sn 1/2j,M 1由邊界條件給出，則用矩陣形式表示為         （26）

20、     第二步，再對(25)式用追趕法求第n 1層的值            摘通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C         本篇論文是由3COME文檔頻道的網(wǎng)友為您在網(wǎng)絡(luò)上收集整理餅投稿至本站的，論文版權(quán)屬原作者，請不用于商業(yè)用途或者抄襲，僅供參考學(xué)習(xí)之用，否者后果自負(fù)，如果此文侵犯

21、您的合法權(quán)益，請聯(lián)系我們。    令Fj=E12ySn 1/2j,l E2LySn 1/2j,l Sn 1/2j,l (1<j<N)    F1=E12ySn 1/21,l E2LySn 1/21,l Sn 1/21,l-C1Sn 10,l    FN=E12ySn 1/2N,l E2LySn 1/2N,l Sn 1/2N,l-C3Sn 1N 1,l    其中Sn 10,l和Sn 1N 1，l由邊界條件給出，則同理可得矩陣

22、方程         (27)這樣，按此步驟一層層地計算。4.3 數(shù)值模擬合理性分析受所掌握的實測資料的限制，目前尚無法對本文提出的算法與含沙量沿垂線分布的實測值進(jìn)行對比。我們用庫里·阿雷克沉沙池的實測資料作了垂線平均值沿程變化的比較。該沉沙池的主數(shù)據(jù)為：池深h=1.53m;平均流速u=0.12m/s;泥沙沉速=0.0176cm/s;懸浮指標(biāo)Z=0.01。計算時取卡門常數(shù)=0.4，a=0.05h。表1給出了計算值和實測值，結(jié)果表明，計算值和實測值比較符合。    表

23、1 斷面平均含沙量驗證(單位：kg/m3)    Verifications of cross-sectional average sediment concentrations            距離(m)    0    200    400    600 &#

24、160;  800    1000    1200            計算值                         

25、       實測值                                    為了進(jìn)一步分析含沙量垂線分布計算結(jié)果的合理性，我們對另一組較粗的泥沙(=0.6

26、16cm/s,Z=0.4)進(jìn)行了對比計算。圖1,圖2分別為兩組沙的計算結(jié)果。從圖中可以看出，計算結(jié)果符合含沙量沿垂線分布的一般規(guī)律，粗沙分布不均勻，細(xì)沙分布較均勻；近底濃度相對較大，水面濃度相對較小，不存在Rouse公式中水面含沙量為0的缺陷；含沙量沿程衰減的特性較為明顯。圖3為較粗一組泥沙的相對含沙量沿垂線分布的沿程變化情況。圖3表明，盡管進(jìn)口斷面按Rouse公式給出了含沙量沿垂線的分布，但由于該斷面實際處于不平衡輸沙狀態(tài)，這種分布并不是穩(wěn)定的。在距進(jìn)口200m處，泥沙的分布調(diào)整到一種不平衡輸沙狀態(tài)，隨著泥沙的沿程淤積，水流輸沙向平衡方向發(fā)展，垂線平均含沙量趨向于水流挾沙力，而含沙量沿垂線分

27、布向平衡時的分布狀態(tài)(Rouse公式)發(fā)展。由于這種發(fā)展是趨向于穩(wěn)定狀態(tài)，因此愈接近下游，分布愈靠近Rouse公式。計算結(jié)果表明，本文提出的計算方法是合理可行的。             圖1 含沙量垂線分布的沿程變化(Z=0.01)            摘通過討論剖面二維非恒定泥沙擴(kuò)散方程的數(shù)值方法，建立了一種用于求解含沙量分布沿程變化的差分格式(Z-C         本篇論文是由3COME文檔頻道的網(wǎng)友為您在網(wǎng)絡(luò)上收集整理餅投