天然河道丁壩群局部沖刷三維數(shù)值模擬

戚 藍(lán)，曾慶達(dá)，吉順文

(1. 天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300072；2. 浙江省水利河口研究院 工程安全研究中心，浙江 杭州 310020)

20世紀(jì)50年代末，灤河開始興建護岸丁壩工程。丁壩具有導(dǎo)流、護岸、防沖和穩(wěn)定河勢的作用，是保灘護岸工程中常見的水工建筑物[1]。但灤河上的丁壩群建設(shè)沒有很好規(guī)劃，丁壩長度、結(jié)構(gòu)形式各異，而且大部分丁壩出現(xiàn)了不同程度的水毀。其中馬良子段以險著稱，若丁壩失效，灤河主流將直沖小埝，危及馬良子村的安全。因此，模擬馬良子丁壩群的局部沖刷，并基于丁壩沖刷規(guī)律對其優(yōu)化十分必要。

目前，丁壩的局部沖刷是眾多專家學(xué)者研究的重點，丁壩建成后由于丁壩上游壅水，會在上游壩根出現(xiàn)一小的立軸回流；至壩頭附近，出現(xiàn)下潛流，產(chǎn)生立軸漩渦并相互碰撞、破碎、向下游運移導(dǎo)致壩頭局部沖刷劇烈。丁壩下游則出現(xiàn)較大的立軸回流及由此于下游壩根處引發(fā)次生反向小回流。這些部位會產(chǎn)生較大的切應(yīng)力，導(dǎo)致河床發(fā)生局部沖刷。喻濤[2]測量了丁壩橫向沖坑的發(fā)展過程。Bey等[3]認(rèn)為丁壩沖坑的產(chǎn)生是由于丁壩的存在使得丁壩附近湍流對河床底部有沖蝕作用，并認(rèn)為最大沖深和沖坑形狀同樣值得研究，因為沖坑形狀一定程度上也會影響丁壩水毀失穩(wěn)。吳學(xué)文等[4]認(rèn)為丁壩局部沖刷主要是壩頭的繞流作用，由此建立了壩頭繞流擠壓流動物理圖示，考慮了非均勻沙的起動問題，得出了非均勻沙河床上丁壩局部沖深公式。事實上，天然河道中的丁壩三維水流情況十分復(fù)雜，許多參數(shù)難以計量，馬永軍[5]分析以往研究結(jié)果發(fā)現(xiàn)，各專家學(xué)者們都只是把回流長度、水流邊界、沖深直接聯(lián)系起來，并沒有找出各參數(shù)之間的關(guān)系。蘇偉等[6]對不同結(jié)構(gòu)形式丁壩進(jìn)行物理模型試驗，發(fā)現(xiàn)丁壩的結(jié)構(gòu)形式對水毀影響明顯。計算機技術(shù)的發(fā)展，為復(fù)雜水流問題的求解提供了計算手段，Tingsanchahi等[7]利用二維平均水深模型并引入校正因子改善k-ε模型，以求解河床底部應(yīng)力。Molls等[8]通過結(jié)合恒定的渦黏性湍流模型改進(jìn)了不穩(wěn)定二維平均水深模型。Jia等[9]利用二維平均水深模型計算丁壩附近水流流線和流速。Duan等[10]用改進(jìn)二維模型對天然河道中泥沙輸運過程進(jìn)行仿真。結(jié)果表明，二維模型只能近似計算泥沙水平運動，無法較好地模擬泥沙在河道中的運輸。越來越多的學(xué)者發(fā)現(xiàn)丁壩附近水流表現(xiàn)出很復(fù)雜的三維特性[11]，因此，開始丁壩的三維數(shù)值模擬研究。Nagata等[11]利用三維水利模型計算了雷諾平均方程和斯托克斯方程，成功模擬了丁壩周圍的泥沙輸運過程。Ouillon等[12]利用三維k-ε模型研究了丁壩周圍三維水流結(jié)構(gòu)和自由液面形狀，求解了壩后沖刷區(qū)水流的水力特性。楊蘭等[13]通過使用Flow-3D軟件，對上挑丁壩群附近流場及局部沖刷進(jìn)行了三維數(shù)值模擬，并得出了可適用于計算丁壩群的湍流模型和推移質(zhì)輸沙率模型。但目前對丁壩群的三維數(shù)值模擬或物理模型模擬大部分針對的是直槽式水槽，而非實際工程中復(fù)雜的天然河道。本文首先基于不同丁壩結(jié)構(gòu)形式的三維水槽模型模擬，對灤河上各種丁壩的沖刷規(guī)律進(jìn)行探索，作為調(diào)整丁壩的依據(jù)。然后對天然河道的丁壩群進(jìn)行三維數(shù)值模擬，對比分析沖刷數(shù)據(jù)與實測地形，其結(jié)果可為工程實際提供一定的參考。

1 理論模型

1.1 控制方程

丁壩群流場三維流動控制方程為連續(xù)性方程和動量方程：

式中：ui為i方向的分速度；u′i為i方向的脈動速度；P為壓力；Sij為應(yīng)變率張量； u′iu′j為雷諾應(yīng)力張量；ρ為流體密度；v為動力黏度；vt為湍流黏度；k為湍動能；δij為克羅內(nèi)克符號(δij=1，i=j，δij=0，i≠j)。

1.2 RNG k-ε方程

標(biāo)準(zhǔn)k-ε模型用于強旋轉(zhuǎn)流或帶有彎曲壁面流動時容易出現(xiàn)失真。RNG k-ε模型通過修正湍動黏性系數(shù)和在ε方程中增加了反映主流的時均應(yīng)變率[7]，可以很好地模擬高應(yīng)變率和流線彎曲程度較大的流動，而天然河道的復(fù)雜地形以及丁壩與水流之間的相互作用會帶來劇烈的水流變形和破碎，且許多專家學(xué)者認(rèn)為RNG k-ε模型對丁壩得到的結(jié)論更為可靠。

RNG k-ε控制方程：

1.3 泥沙輸運模型

目前，研究丁壩局部沖刷的規(guī)律主要有5個理論途徑，分別是流速、切應(yīng)力、能量平衡、統(tǒng)計法則及沙波運動，而湍流流場和泥沙輸運對于丁壩群局部沖刷坑的發(fā)育影響在理論層面還不是十分成熟?；谇叭藢Χ螞_刷規(guī)律的試驗以及泥沙運動理論分析，選用以希爾茲(Shields)數(shù)為基礎(chǔ)的泥沙推移質(zhì)輸沙率模型作為此次丁壩群沖刷三維數(shù)值模擬的理論模型。

與楊蘭等[13]的泥沙輸運模型相同，模型計算中不考慮懸移質(zhì)對沖刷坑發(fā)育的影響，只分析河床底部推移質(zhì)的沖刷作用，并認(rèn)為推移質(zhì)為均勻的泥沙顆粒。采用的推移質(zhì)輸沙率公式以切應(yīng)力為主要參數(shù)，且切應(yīng)力與輸沙率的關(guān)系成正比。

水流切應(yīng)力:

泥沙起動的臨界切應(yīng)力:

希爾茲數(shù):

臨界希爾茲數(shù):

式中：U*為摩阻流速；U*c為臨界摩阻流速；ρ和ρs分別為水流和泥沙密度；d為泥沙平均粒徑。假定單寬推移質(zhì)輸沙率計算式為：

式中：qb為單寬推移質(zhì)輸沙率；ub為推移質(zhì)的平均輸運速度；p為泥沙起動概率。對于實際工程地質(zhì)情況的土工試驗，可以根據(jù)上式得到ub和p，即可求解qb。

動摩擦力為:

FD=fD并聯(lián)立(7)，得到:

式中：aU*為推移質(zhì)運動的水流速度，當(dāng)離沙床較近時，a=6～10；CD為推力系數(shù)。當(dāng)θ0=0時，ub=0，θ0相當(dāng)于止動相對切應(yīng)力，應(yīng)小于臨界相對切應(yīng)力θc。由此，式(15)可寫為:

取泥沙臨界切應(yīng)力和運動顆粒所受切應(yīng)力之和為切應(yīng)力，即:

式中：n為泥沙顆?？倲?shù)，而單位面積河床泥沙顆?？倲?shù)1/d2與p的關(guān)系為：p=n/d-2。聯(lián)立式(9)和(10)，得到床面泥沙顆粒的起動概率為

綜上所述，影響單寬推移質(zhì)輸沙率的兩個重要因素為ub和p，將ub和p的表達(dá)式代入(11)，并取動摩擦系數(shù)β=0.08，常數(shù)a=9.3?？傻萌缦峦埔瀑|(zhì)輸沙公式：

式(19)即為所采用的泥沙運動模型?？梢园l(fā)現(xiàn)，對于某些泥沙顆粒，推移質(zhì)輸沙率只與希爾茲數(shù)有關(guān)，只有當(dāng)希爾茲數(shù)大于臨界希爾茲數(shù)時，泥沙才會起動被沖刷，采用θc=0.05。

2 模型驗證

實際工程丁壩段如圖1所示，可見11#丁壩相比周圍的丁壩明顯較長，不夠合理。這是由于丁壩有束窄河床的作用，河床被束窄后，河道斷面面積減小；雖然丁壩前有一定程度的壅水，但斷面面積仍有大比例收縮，導(dǎo)致相同流量下流速迅速增大，特別是河道底部的切應(yīng)力會迅速增大，下切河床造成沖刷。同時，該河段自9#丁壩后河道寬度明顯縮窄，流速更大，由于灤河段河床地質(zhì)條件差，基本是粉細(xì)砂組成，更容易遭到?jīng)_刷。

對比河底實測地形圖發(fā)現(xiàn)，11#丁壩后有很深的沖刷坑，最深處超過13 m，因此選擇11#丁壩作為主要的研究對象，研究不同丁壩長度對丁壩群的沖刷影響。

計算模型長1 400 m，寬1 000 m，高30 m，其中11#丁壩上游400 m至下游300 m為動床計算區(qū)域，其余為定床。河床主槽糙率為0.02，平均水深為5 m，為非淹沒丁壩群，計算丁壩群位于左岸，計算河道模型平面輪廓與天然河道一致，考慮到河道寬度遠(yuǎn)遠(yuǎn)大于水深，因此河道斷面形狀取矩形。原始11#丁壩長為34.5 m，動床所鋪泥沙厚度為20 m，泥沙平均粒徑d=0.018 mm，密度ρs=2 650 kg/m3。計算區(qū)域采用矩形網(wǎng)格劃分為A，B，C區(qū)，A區(qū)和C區(qū)為定床，B區(qū)為動床。對動床部分加密網(wǎng)格，總網(wǎng)格數(shù)量約133萬，網(wǎng)格分塊劃分見圖2。

計算時采用有限體積法離散控制方程，湍流模型為RNG k-ε模型。上游設(shè)為速度進(jìn)口，根據(jù)上游實測流速和平均深度，給定模型流量為8 430 m3/s，平均水深5 m，出口設(shè)為自由出流，上表面采用VOF法捕捉液體自由面，固體邊界為無滑移面，近壁面采用標(biāo)準(zhǔn)壁面函數(shù)進(jìn)行處理。對比11#壩頭附近的沖刷深度與實測地形圖，實測與模擬沖刷地形擬合較好，如圖3所示?？梢?，在復(fù)雜天然河道形狀下，可以利用RNG k-ε模型和推移質(zhì)輸沙率模型模擬丁壩的局部沖刷。

圖 1 馬良子段丁壩群布置Fig. 1 Layout of spur dike group along Maliangzi reach

圖 2 計算區(qū)域網(wǎng)格剖分Fig. 2 Meshes of computational domain

圖 3 壩頭附近沖刷深度橫縱斷面對比Fig. 3 Comparison of transverse and longitudinal sections of the scouring depth near the spur dike head

3 基于丁壩壩長的調(diào)整方案

灤河上馬良子段丁壩結(jié)構(gòu)形式簡單，多為下挑式丁壩，挑角約70°～85°，變化幅度不大，且壩頭形式單一，為直立式丁壩。通過丁壩設(shè)計資料以及地形勘察資料發(fā)現(xiàn)，馬良子段11#丁壩設(shè)計明顯過長，設(shè)計長度34.5 m，比相鄰丁壩長約10 m，而且10#～12#丁壩在河流收縮段。測量數(shù)據(jù)也表明11#丁壩壩后是沖刷最嚴(yán)重的區(qū)域，最大沖深達(dá)13.98 m。因此丁壩長度不合理很可能是造成嚴(yán)重沖刷的原因。本文將原壩長方案定為方案1，將11#丁壩縮短5 m和10 m，對應(yīng)方案2和3，對丁壩周圍水力特性及沖刷特性進(jìn)行研究。

3.1 丁壩繞流流場與流速分布

方案1～3計算域水流流場與流速如圖4所示。可見，隨著11#丁壩壩長的減小，水流隨之平順，丁壩后端壩田附近形成的漩渦減小。同時，在原有丁壩群間距不變的情況下，縮短11#丁壩依然可以保持丁壩群對水流的有效影響，即壩前水流流入壩后壩田區(qū)與壩田內(nèi)水流相互作用，與岸邊不直接接觸，主河槽流速較兩岸大，方案1中過11#壩后主河槽流速明顯增加，方案2和3增速不明顯。

圖 4 11#丁壩附近繞流流場及流速分布Fig. 4 Flow field and velocity distribution around No.11 spur dike

3.2 沖刷分析

丁壩群沖刷形態(tài)對比如圖5所示，可以發(fā)現(xiàn)調(diào)整方案1和2中丁壩群最大沖深坑出現(xiàn)在11#丁壩壩頭附近，并且在主流區(qū)形成一狹長沖刷帶；而方案3丁壩群則沒有明顯沖刷坑，沖刷深度較為一致。這說明縮短丁壩群中過長的丁壩可以有效減輕其壩頭的局部沖刷，丁壩群調(diào)整后河床沖刷明顯減輕。

圖 5 丁壩群沖刷形態(tài)對比Fig. 5 Comparison of scour patterns around spur dikes

根據(jù)實測沖刷地形在11#丁壩壩頭附近設(shè)置測點，其中測點1是原壩后沖坑最深處，測點2和4分別是測點1上、下游100 m位置，測點3為測點1的右邊5 m處，測點布置如圖6所示。丁壩群沖刷數(shù)值對比如表1所示?？梢?，縮短過長的丁壩長度可以有效減小丁壩局部沖刷深度，從而降低丁壩的水毀風(fēng)險。

圖 6 測點布置Fig. 6 Layout of measuring points

表 1 3種方案各測點數(shù)值對比Tab. 1 Comparison of values of measuring points in three schemes

4 結(jié) 語

基于天然河道形狀，利用推移質(zhì)輸沙率模型和VOF自由表面處理法，對非淹沒丁壩群的繞流流場、流速、局部沖刷進(jìn)行三維數(shù)值模擬研究，并與實測資料進(jìn)行對比，得到以下主要結(jié)論：丁壩群中某一丁壩長度較周圍丁壩明顯偏長易引起河道水流集中，壩后漩渦加劇，且丁壩長度愈長，阻水能力愈大，壩頭沖刷越劇烈。最大沖刷發(fā)生在壩頭處，且在沖坑內(nèi)略偏壩軸線下游，上游部分坡度較陡，范圍較小，下游部分坡度較緩，范圍較大。說明丁壩群設(shè)計時應(yīng)合理規(guī)劃丁壩長度，丁壩群內(nèi)各個丁壩長度不協(xié)調(diào)可能會形成新的水害。

由于實際工程中丁壩的沖刷坑形成非常復(fù)雜，本次模擬的沖刷坑下游沖刷與實際還有一些偏差，對丁壩局部沖刷機理和實際水流過程有待進(jìn)一步研究。