內(nèi)河航道挖槽立面二維水沙數(shù)學(xué)模型研究

吳 騰，吳玲莉，丁 飛

（河海大學(xué)港口海岸與近海工程學(xué)院，南京210098）

我國(guó)河流眾多，水運(yùn)資源十分豐富，但大部分河流含沙量較大，航道泥沙淤積嚴(yán)重，不得不消耗大量資金進(jìn)行挖槽清淤以改善航道條件。由于航道挖槽改變了原有的河床形態(tài)，使水流流態(tài)更加復(fù)雜，同時(shí)也破壞了原有的水沙間的相對(duì)平衡，如何確定挖槽斷面的尺寸，使挖槽后的泥沙回淤量較小就成為航道挖槽的關(guān)鍵問(wèn)題。

關(guān)于減小航道淤積的研究由來(lái)以久，但由于航道疏浚問(wèn)題的復(fù)雜性，且相關(guān)實(shí)測(cè)資料較少，使得該問(wèn)題的研究難度較大。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，許多學(xué)者也采用數(shù)學(xué)模型研究該問(wèn)題［1］。黃永健結(jié)合沿水流方向的一維水流運(yùn)動(dòng)和垂向二維泥沙擴(kuò)散方程求解斷面平均流速和含沙量濃度分布，再進(jìn)行水流穿越河槽的不平衡輸沙計(jì)算，得到航道挖槽的回淤量［2］，由于在求解過(guò)程中采用了一維模型計(jì)算流速，難以反映挖槽后由于河床變形產(chǎn)生的環(huán)流對(duì)泥沙輸移的影響，因此，在不平衡輸沙計(jì)算中挾沙力的計(jì)算結(jié)果與理論存在一定差異；為了更好研究航道挖槽后泥沙的回淤情況，夏軍強(qiáng)建立剖面二維水流泥沙數(shù)學(xué)模型，該模型采用流函數(shù)-渦量方程，可以較好克服N-S 方程中計(jì)算動(dòng)水壓強(qiáng)較難的問(wèn)題［3］。對(duì)于挖槽立面二維水沙數(shù)學(xué)模型而言，自由面和河床形態(tài)是在不斷變化的，隨著計(jì)算時(shí)間的推移，模型的計(jì)算區(qū)域都可能隨時(shí)間變化，如果保持網(wǎng)格不變化，在計(jì)算過(guò)程中難免帶來(lái)誤差，當(dāng)誤差大到一定程度就有可能導(dǎo)致程序的振蕩與發(fā)散，對(duì)模型計(jì)算極為不利。因此，在模型中應(yīng)充分考慮自由水面和河床變化問(wèn)題。

水體產(chǎn)生的壓強(qiáng)常簡(jiǎn)化為靜水壓強(qiáng)，在航道挖槽后，河床地形發(fā)生劇烈變化，水流流態(tài)也發(fā)生較大變化，對(duì)于立面二維模型，需要考慮動(dòng)水壓強(qiáng)與靜水壓強(qiáng)的差異。為此，本文擬采用有限體積法建立內(nèi)河航道挖槽回淤立面二維水沙數(shù)學(xué)模型。該模型將采用動(dòng)網(wǎng)格技術(shù)，使計(jì)算區(qū)域隨流動(dòng)區(qū)域的變化而改變，這將有效地減少網(wǎng)格與計(jì)算區(qū)域不重合帶來(lái)的誤差，更好地模擬自由表面和河床的變化；同時(shí)，模型計(jì)算中將引入靜水壓強(qiáng)修正值，克服因挖槽斷面形態(tài)變化引起的流速劇烈變化帶來(lái)的模擬困難。采用該模型探討航道挖槽后泥沙回淤的機(jī)理，分析了不同挖槽斷面形態(tài)與流態(tài)變化和泥沙淤積的關(guān)系，為挖槽尺寸的優(yōu)化設(shè)計(jì)提供參考。

1 內(nèi)河航道挖槽立面二維模型的建立

1.1 基本方程

航道挖槽后，槽內(nèi)水沙運(yùn)動(dòng)可采用下列方程進(jìn)行描述［4］。

水流連續(xù)方程

水流運(yùn)動(dòng)方程

懸移質(zhì)運(yùn)動(dòng)方程

河床變形方程

式中：u 為沿水流方向流速；w 為垂向流速，以向上方向?yàn)檎籶 為動(dòng)水壓強(qiáng)；ρ 為清水密度；vt為紊流粘滯性系數(shù)；h 為水深；ε 為泥沙紊動(dòng)擴(kuò)散系數(shù)；Zb為河床高程；ω 為泥沙沉速；pr為計(jì)算參數(shù)；sa、s*a床面含沙量及挾沙力；γ′泥沙干容重。

1.2 考慮動(dòng)水壓強(qiáng)的水流方程

（1）動(dòng)水壓強(qiáng)的引入。

航道進(jìn)行挖槽后，河道斷面形態(tài)發(fā)生較大改變，水體的動(dòng)水壓強(qiáng)分布不能簡(jiǎn)化為靜水壓槍分布，為了使水位函數(shù)ξ（x，t）與動(dòng)量方程密切聯(lián)系起來(lái)，將動(dòng)水壓強(qiáng)分解為靜水壓強(qiáng)和一個(gè)壓強(qiáng)修正值

式中：p′為壓強(qiáng)修正值，是在水體流動(dòng)時(shí)由于流線(xiàn)彎曲和流速不均勻所產(chǎn)生的附加壓強(qiáng)；ps為靜水壓強(qiáng)，可表示為

式中：ξ 為水位；y 為河床高程；ρ 為清水密度；g 為重力加速度。

將壓強(qiáng)分解后，壓強(qiáng)梯度可分別改為

將式（9）代入水流運(yùn)動(dòng)方程式（2）和式（3），水流運(yùn)動(dòng)方程變?yōu)?/p>

（2）自由表面的處理。

式中：u 為x 方向沿水深平均流速；h 為水深；ξ 為水位；vs為水面流速。在程序編制中，可事先預(yù)留部分網(wǎng)格，判斷是否過(guò)水來(lái)確定是否進(jìn)入程序計(jì)算。

（3）方程統(tǒng)一形式。

綜合上述修改，方程（1）、（10）、（11）、（4）、（5）和（12）即為本文建立的航道挖槽回淤立面二維水沙數(shù)學(xué)模型基本方程，連續(xù)方程和動(dòng)量方程可寫(xiě)成統(tǒng)一形式

方程的離散采用有限體積法，離散形式可參考文獻(xiàn)［5］。

表1 式（13）各參變量形式Tab.1 Parameters expression in formula（13）

2 相關(guān)問(wèn)題的處理

（1）進(jìn)口斷面。

進(jìn)口含沙量分布依據(jù)張瑞瑾、丁君松方法的含沙量垂向分布［6］

式中：s 為進(jìn)口斷面平均含沙量；ξ 為相對(duì)水深；I 為參數(shù)，根據(jù)丁君松的研究，I 與懸浮指標(biāo)Z 的關(guān)系應(yīng)為

（2）水面泥沙邊界條件。

（3）河底泥沙邊界條件。

式中：Sb*為近底挾沙力，文中采用van Rijn 提出的方法進(jìn)行計(jì)算［7］。

（4）關(guān)于運(yùn)動(dòng)粘性系數(shù)vt的確定。

文中紊動(dòng)粘滯性系數(shù)由下式確定［8］

式中：U*為摩阻流速；H 為斷面平均水深；α 為常數(shù)，α=0.25～1.0。

3 模型的驗(yàn)證及分析

3.1 流速的驗(yàn)證

模型的驗(yàn)證采用文獻(xiàn)［3］中的挖槽試驗(yàn)數(shù)據(jù)，試驗(yàn)的尺寸及水流條件為槽溝上下邊坡坡度均為1:2，槽底寬為1 m，深度為0.2 m。進(jìn)口處水深為0.2 m，垂線(xiàn)平均流速為0.4 m/s。同時(shí)假定在進(jìn)口處于輸沙平衡狀態(tài),即床面保持不沖不淤，相應(yīng)的懸移質(zhì)垂線(xiàn)平均含沙量為0.12 kg/m3，中值粒徑為0.1 mm。計(jì)算網(wǎng)格為50×20個(gè)網(wǎng)格，網(wǎng)格長(zhǎng)度分別為0.05 m 和0.025 m。圖1 為0.4 m、1.1 m、2.2 m 處斷面垂向流速驗(yàn)證。

可以看出進(jìn)口斷面與出口斷面流場(chǎng)較為平順，在挖槽由于斷面形態(tài)突變，流場(chǎng)較為復(fù)雜。定量上，模型的計(jì)算值與實(shí)測(cè)值較為接近，能反映挖槽斷面的流速變化規(guī)律，可用于挖槽的研究。圖2 為挖槽斷面流場(chǎng)圖，圖3 為挖槽坡腳局部放大流場(chǎng)圖，受挖槽斷面形態(tài)影響，挖槽坡腳出現(xiàn)明顯的回流。

3.2 垂線(xiàn)含沙量計(jì)算分析

圖4 為挖槽不同斷面垂線(xiàn)含沙量分布，含沙量濃度分布上稀下濃，進(jìn)口斷面含沙量上下濃度差異較大，隨著水深的增大，挖槽斷面逐漸擴(kuò)大，產(chǎn)生回流，斷面上下濃度交換增多，斷面的含沙量梯度減小，挖槽后的斷面流速逐漸平順，與河床發(fā)生交換，垂向濃度又稍有增大，與定性分析相同，說(shuō)明該模型能計(jì)算含沙量的分布。

4 結(jié)論

本文將模型中的壓強(qiáng)項(xiàng)分解為靜水壓強(qiáng)和修正壓強(qiáng)項(xiàng)，建立了動(dòng)水壓強(qiáng)的立面二維水沙數(shù)學(xué)模型，同時(shí)該模型采用了動(dòng)網(wǎng)格技術(shù)，使計(jì)算區(qū)域隨流動(dòng)區(qū)域的變化而改變。采用試驗(yàn)資料對(duì)該模型進(jìn)行檢驗(yàn)，結(jié)果表明，該模型能較好反應(yīng)航道挖槽后流速、含沙量的變化，能清晰模擬由挖槽產(chǎn)生的橫軸環(huán)流；此外，模型還能較好模擬航道挖槽后斷面含沙量的垂向分布，可供挖槽后的回淤估算提供參考。

［1］Leedertes J J. A Water Quality Simulation Model for Well-Mixed Estuaries and Coastal Seas［M］. Santa Monica:Rand，1970.

［2］黃永健. 長(zhǎng)江口挖槽自然回淤的計(jì)算［J］. 泥沙研究，1997（2）：69-73.HUANG Y J. Calculation of Natural siltation in Yangtze Delta excavation［J］.Sediment Research，1997（2）：69-73.

［3］夏軍強(qiáng)，談廣鳴. 橫向槽溝內(nèi)泥沙淤積與水平軸環(huán)流變化的數(shù)值模擬［J］.水利學(xué)報(bào)，1998（8）:51-56.XIA J Q，TAN G M. Numerical Simulation of sediment deposition and horizontal axis circulation flow in lateral trenches［J］.Journal of hydraulic engineering，1998（8）:51-56.

［4］余明輝，吳騰，楊國(guó)錄. 剖面二維水沙數(shù)學(xué)模型及其初步應(yīng)用［J］.水力發(fā)電學(xué)報(bào)，2006，25（4）:66-69.YU M H，WU T，YANG G L.Study on vertical 2-D sediment numerical model and its primary application［J］.Journal of hydroelectric engineering，2006，25（4）:66-69.

［5］吳騰. 壩區(qū)水沙立面二維數(shù)學(xué)模型研究［D］. 武漢：武漢大學(xué)，2005.

［6］張瑞瑾. 河流泥沙動(dòng)力學(xué)［M］.北京：中國(guó)水利電力出版社，1998.

［7］Van Rijn L C. Sediment transport，part Ⅱ：suspended load transport［J］. Journal of hydraulic Engineering，ASCE，1984，110（11）:1 613-1 641.

［8］韓其為. 黃河泥沙若干理論問(wèn)題研究［M］.鄭州：黃河水利出版社，2010.