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

1、剖面二維非恒定懸移質(zhì)泥沙擴(kuò)散方程的數(shù)值方法1 引言 數(shù)學(xué)模擬方法正在成為 研究 河流泥沙 問題 的重要手段。 目前 ，一維數(shù)學(xué)模型 發(fā)展 較成熟，已廣泛 應(yīng)用 于模擬長河段的長期變形，但它只能給出河段平均沖淤深度的沿程變化，如需了解短河段的河床變形細(xì)節(jié)，則要采用二維以至三維數(shù)學(xué)模型。不論是一維數(shù)學(xué)模型還是平面二維維數(shù)學(xué)模型，都不能反映含沙量沿垂線的分布狀況，并忽略了含沙量沿垂線分布對垂線平均含沙量變化過程的 影響 。要解決這類問題，必須建立剖面二維數(shù)學(xué)模型。這種模型主要通過解剖面二維泥沙擴(kuò)散方程來研究懸移質(zhì)泥沙沿水深的分布及含沙量的變化過程，對水電站進(jìn)口和其它引水工程的引水口高程的確定都能提供

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

3、 為水流方向和鉛直方向的維軸； u,v 分別為沿水流方向和鉛直方向的時(shí)均流速； sx , sy 分別為水流方向和鉛直方向的泥沙擴(kuò)散系數(shù)； ,S 分別為泥沙靜水沉速和含沙量。 對于式 (1) 的求解，研究者一般會對它進(jìn)行不同程度的簡化，為此引入以下假定中的一種或幾種： A. 非恒定流可以概化為梯級式恒定流，即 ; B. 在一個(gè)時(shí)段內(nèi)，認(rèn)為泥沙運(yùn)動可以概化為處于恒定狀態(tài)，即 ; C. 在二維流動中，縱向擴(kuò)散系數(shù)與方程其他項(xiàng)相比，可以忽略不計(jì)，即認(rèn)為方程右端第一項(xiàng)可以忽略； D. 認(rèn)為懸移質(zhì)泥沙粒徑均一，即 =const; E. 認(rèn)為水流為二維均勻流，即 v=0 。 為簡單起見，我們討論的范圍限于水

4、流條件為二維非恒定均勻流，懸移質(zhì)泥沙粒徑均勻，為此引入假定 C 、 D 、 E 。這時(shí)，泥沙擴(kuò)散方程為 (2) 目前，對 s 的變化 規(guī)律 研究得不很充分，一般假定 s = m (3) 其中 m 為動量傳遞系數(shù)，為修正值。由勃蘭特爾摻長 理論 可得 m = u*y/(1-y/h) (4) 式中為卡門常數(shù)， u* 為摩阻流速。 對于 u ，我們?nèi)】?- 勃蘭特爾對數(shù)流速分布公式 (umax-u)/u*=(1/ )ln(h/y) (5) 令 W= + ( u*/h)(1-2y/h), 則式 (2) 可變形為 (6) 3 差分方程 3.1 網(wǎng)格的剖分 為建立差分方程，首先必須剖分網(wǎng)格。我們?nèi)r(shí)間步

5、長 t= ， X 方向的空間步長 x=h1 ， Y 方向的空間步長為 y=h2, 這樣形成如下網(wǎng)格 Dh=(xj,yl,tn) xj=jh1,yl=lh2,tn=n 其中 j=0,.,N; l=0,.,M;n 0; 3.2 構(gòu)造差分格式 通過對流方程和擴(kuò)散方程的差分格式的構(gòu)造，我們可以得到對流擴(kuò)散方程的差分格式。由于隱式格式穩(wěn)定性好，考慮 Crank-Nicholson 型隱式格式。為此，引入差分算子記號 2 y Snj,l=Snj,l+1-2Snj,l+Snj,l-1 LxSnj,l=Snj+1,l-Snj-1,l LySnj,l=Snj,l+1-Snj,l-1 為了看得更清楚，暫且取 h1

6、=h2=h. 對式 (6) 離散，則 C-N 格式為 Sn+1j,l-Snj,l/ +u/4hLx(Snj,l+Sn+1j,l)= s /2h2 2 y (Sn+1j,l+Snj,l)+W/4hLy(Sn+1j,l+Snj,l) (7) C-N 格式的精度是二階的，絕對穩(wěn)定。但對于二維 問題 ，由 (7) 導(dǎo)出的方程組，其系數(shù)矩陣不是三對角矩陣，不能用追趕法求解。因此，考慮構(gòu)造交替方向的隱式格式 ( 命名為 Z-C 格式 ) (Sn+1/2j,l-Snj,l)/( /2)+u/2hLxSnj,l= s /h2 2 y S n+1/2 j,l +W/2hLySn+1/2 j,l (8) (Sn+

7、1j,l-Sn+1/2j,l)/ /2+u/2hLxSn+1j,l= s /h2 2 y Sn+1/2j,l+W/2hLySn+1/2j,l (9) 可以看出， 計(jì)算 Sn+1j,l 是由兩步組成的，每一步僅是一個(gè)方向的隱式，故用兩次追趕法即可。 3.3 精度 分析 現(xiàn)在，我們考慮 Z-C 格式的精度。先設(shè)法消去過渡值 Sn+1/2j,l ，為此，將 (8) 和 (9) 兩式相加，可得 Sn+1j,l-Snj,l+ u/2hLx(Sn+1j,l+Snj,l)= s /2h2 2 y Sn+1/2j,l+ W/4hLxSn+1/2j,l (10) 將 (8) 和 (9) 兩式相減，可得 Sn+1

8、/2j,l=(Sn+1j,l+Snj,l)/2+ s /4h2 2 y (Snj,l-Sn+1j,l)+ W/8hLy(Snj,l-Sn+1j,l) (11) 把式 (11) 代入 (10) ，變形整理，可得 (Sn+1j,l-Snj,l)(1- u s /8h3Lx 2 y - 2 uW/16h2L x Ly)=(Snj,l+Sn+1j,l)( u s /2h2 2 y + 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,

9、tn+1)-S(xj,yl,tn)/ (1- 2 s /8h3Lx 2 y - 2 uW/16h2LxLy)-S(xj,yl,tn+1)-S(xj,yl,tn)/ ( s /2h2) 2 y + W/4hLy- u/4hLx)=O( 2 +h2) (13) 由此可見， Z-C 格式具有二階精度。 3.4 穩(wěn)定性分析 現(xiàn)在，我們來討論 Z-C 格式的穩(wěn)定性。為此，把式 (12) 變形整理得 (1- s /2h 2 2 y - W/4hLy)(1+ u/4hLx)Sn+1j,l=(1+ s /2h 2 y 2 + W/4hLy)(1- u/4hLx)S n j,l (14) 由式 (14) 可得出

10、過渡因子為 G( ,k)=(1-2 s /h 2 sin2k2h/2+i W/2hsink2h)(1-i u/2hsink1h)/(1+2 s /h2sin2k2h/2-i W/2hsink2h)(1+i u/2hsink1h) (15) 令 a=2 s /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)|

11、21, 所以 Z-C 格式是絕對穩(wěn)定的。 我們考慮初邊值 問題 。 (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.01 0.05h 。 (2) 水面條件 (y=h) (18) (3) 底部邊界條件 (y=a) (19) 式中 Sa* 為近底挾沙力，即輸沙平衡時(shí)的近底售沙量 Sa. (4) 近底含沙量計(jì)算近底含沙量在求解泥沙擴(kuò)散方程時(shí)具有邊界條件性質(zhì)，它選取的正確與否，意味著所給邊界條件是否正確。實(shí)際工程中一般缺乏實(shí)測資料，近底

12、含沙量不易測定。這里，我們利用水流挾沙力和含沙量沿垂線分布公式來反求近底含沙量 。 已知斷面平均挾沙力為 S*=k(u3/gh )m (20) 假定 Sa*= S* (21) 輸沙平衡時(shí)，含沙量沿垂線分布用 Rouse 公式 (17) 表示，用 (17) 表達(dá)挾沙力的垂線分布 ，然后沿垂線積分得斷面平均挾沙力為 (22) 將 (22) 與 (21) 比較，可得 (23) 4.2 計(jì)算步驟 為方便計(jì)算，將式 (8) 和 (9) 式變形整理，并對 X ， Y 方向取不同的空間步長 ( s /2h22- W/4h2)Sn+1/2j,l+1-( s /h22+1)Sn+1/2j,l+( s /2h22

13、- W/4h2)Sn+1/2j,l+1= u/4h1LxSnj,l-Snj,l (24) - u/4h1Snj-1,l+Snj,l+ u/4h1Snj+1,l= s /2h22 2 y Sn+1/2j,l+ W/4h2LySn+1/2j,l+Sn+1/2j,l (25) 在一個(gè)時(shí)間層 ( 第 n 層 ) 內(nèi)，計(jì)算分兩步進(jìn)行： 第一步，對式 (24) 用追趕法求第 n+1/2 層的過渡值。 令 C1=- u/4h1,C2=1,C3= u/4h1,E1= s /2h22,E2= W/4h2 D1= s /2h22- W/4h2,D2=- s /h22-1,D3= s /2h22+ W/4h2 l=

14、1 時(shí)， D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1 l=2 時(shí)， D1Sn+1/2j,0+D2Sn+1/2j,1+D3Sn+1/2j,2=C3LxSnj,1-Snj,1 l=M 時(shí) ,D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3L x Snj,M-Snj,M 令 Hl=-Snj,l+C3LxSnj,l (1lM) 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+

15、1 由邊界條件給出，則用矩陣形式表示為 （ 26 ） 第二步，再對 (25) 式用追趕法求第 n+1 層的值 令 Fj=E1 2 y Sn+1/2 j,l +E2LySn+1/2j,l+Sn+1/2j,l (1jN) F1=E1 2 y Sn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C 1 S n+1 0,l FN=E1 2 y Sn+1/2N,l+E2L y Sn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l 其中 Sn+10,l 和 Sn+1N+1 ， l 由邊界條件給出，則同理可得矩陣方程 (27) 這樣，按此步驟一層層地計(jì)算。 4.3 數(shù)值模擬合理性

16、分析 受所掌握的實(shí)測資料的限制， 目前 尚無法對本文提出的算法與含沙量沿垂線分布的實(shí)測值進(jìn)行對比。我們用庫里阿雷克沉沙池 的實(shí)測資料作了垂線平均值沿程變化的比較。該沉沙池的主要數(shù)據(jù)為：池深 h=1.53m; 平均流速 u=0.12m/s; 泥沙沉速 =0.0176cm/s; 懸浮指標(biāo) Z=0.01 。計(jì)算時(shí)取卡門常數(shù) =0.4 ， a=0.05h 。表 1 給出了計(jì)算值和實(shí)測值，結(jié)果表明，計(jì)算值和實(shí)測值比較符合。 表 1 斷面平均含沙量驗(yàn)證 ( 單位： kg/m3) Verifications of cross-sectional average sediment concentrations

17、 距離 (m) 0 200 400 600 800 1000 1200 計(jì)算值 3.012 2.573 2.071 1.639 1.209 0.849 0.539 實(shí)測值 3.012 2.650 1.810 1.520 1.141 0.822 0.561 為了進(jìn)一步分析含沙量垂線分布計(jì)算結(jié)果的合理性，我們對另一組較粗的泥沙 ( =0.616cm/s,Z=0.4) 進(jìn)行了對比計(jì)算。圖 1, 圖 2 分別為兩組沙的計(jì)算結(jié)果。從圖中可以看出，計(jì)算結(jié)果符合含沙量沿垂線分布的一般 規(guī)律 ，粗沙分布不均勻，細(xì)沙分布較均勻；近底濃度相對較大，水面濃度相對較小，不存在 Rouse 公式中水面含沙量為 0 的缺

18、陷；含沙量沿程衰減的特性較為明顯。圖 3 為較粗一組泥沙的相對含沙量沿垂線分布的沿程變化情況。圖 3 表明，盡管進(jìn)口斷面按 Rouse 公式給出了含沙量沿垂線的分布，但由于該斷面實(shí)際處于不平衡輸沙狀態(tài)，這種分布并不是穩(wěn)定的。在距進(jìn)口 200m 處，泥沙的分布調(diào)整到一種不平衡輸沙狀態(tài)，隨著泥沙的沿程淤積，水流輸沙向平衡方向 發(fā)展 ，垂線平均含沙量趨向于水流挾沙力，而含沙量沿垂線分布向平衡時(shí)的分布狀態(tài) (Rouse 公式 ) 發(fā)展。由于這種發(fā)展是趨向于穩(wěn)定狀態(tài)，因此愈接近下游，分布愈靠近 Rouse 公式。計(jì)算結(jié)果表明，本文提出的計(jì)算 方法 是合理可行的。 圖 1 含沙量垂線分布的沿程變化 (Z=0.01) 圖 2 含沙量垂線分布的沿程變化 (Z=0.4) Changes of vertical distributions of sediment concentrations(Z=0.01) Chang