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

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

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

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

1 理論模型

1.1 控制方程

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

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

1.2 RNG k-ε方程

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

RNG k-ε控制方程：

1.3 泥沙輸運(yùn)模型

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

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

水流切應(yīng)力:

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

希爾茲數(shù):

臨界希爾茲數(shù):

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

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

動(dòng)摩擦力為:

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

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

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

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

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

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

2 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

3.1 丁壩繞流流場(chǎng)與流速分布

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

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

3.2 沖刷分析

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

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

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

圖 6 測(cè)點(diǎn)布置Fig. 6 Layout of measuring points

表 1 3種方案各測(cè)點(diǎn)數(shù)值對(duì)比Tab. 1 Comparison of values of measuring points in three schemes

4 結(jié) 語(yǔ)

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

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