劉玉嬌 余明輝 田浩永

(武漢大學水資源與水電工程科學國家重點實驗室,武漢 430072)

引言

天然河道的演變是水流攜帶泥沙并與河道邊界相互交換相互影響的辯證統(tǒng)一體.其中,縱向流速分布規(guī)律與河道過流能力[1]和泥沙輸移[2]的密切相關(guān),也與崩岸、切灘等突發(fā)現(xiàn)象相關(guān).對于河流流速的分布規(guī)律,研究主要從模型試驗[3-5]、數(shù)值模擬[6-9]和解析解計算[10-13]進行.解析解由于效率高,方法簡單,在實際應用方面被廣泛采用[14].在順直復式河道中,Shiono 和Knight[15]提出了SKM 方法,認為二次流導致的切應力沿橫斷面為線性分布,對動量方程沿水深積分,推導出垂線平均縱向流速(以下簡稱為Ud,下標d 表示為垂線平均量)的解析解;Devi 等[16]進一步根據(jù)實測數(shù)據(jù)量化二次流中的參數(shù),發(fā)現(xiàn)參數(shù)與定水深和變化的水深之間有不同的關(guān)系;復式河槽漫灘區(qū)域常為順直的,但主河槽則傾向于彎曲的,存在彎道環(huán)流的影響,Ervine 等[17]考慮了復式河道彎曲主槽中水流的強三維性,認為二次流項與Ud的二次方成正比關(guān)系,對SKM 方法進行了改進;Knight等[18]分析了環(huán)流結(jié)構(gòu)與二次流項取值正負之間的關(guān)系,認為環(huán)流中心左側(cè)和右側(cè)二次流項符號相反;Liu 等[19]和Huai 等[20]則進一步在動量方程中增加了灘地植被影響,在解析解公式中考慮了植被對流速的影響;Zhong 等[21]考慮了在兩岸及水面結(jié)冰條件下,解析解的形式,并分析了床面阻力、冰蓋阻力和二次流等因素的影響;楊中華和高偉[22]比較了三種不同的積分假定,將計算值與實測值比較,得出最為吻合的計算公式;許唯臨[11]將二次流項歸并入雷諾切應力項,對復式順直河道和復式游蕩型河道中漫灘流速進行了研究,復式游蕩型河道流速計算其采用了分層計算的方法.哈岸英等[23]利用SKM 方法研究了沖積河流漫灘水流水力因子隨水深的變化特性,提出了定量劃分復式斷面灘與槽的判別指標.Rezaei 和Knight[24]探究了SKM 方法在非對稱復合水槽中的應用.Sharifi 和Sterling[2]利用SKM 方法計算床面切應力的橫向分布,并將其與泥沙輸移公式進行的對比,以獲得對泥沙輸移的預測.

彎曲型河道是常見的天然河道形態(tài),水流經(jīng)過彎道時受到離心力的作用,產(chǎn)生橫向壓力梯度,形成彎道橫向環(huán)流,橫向環(huán)流與縱向流速結(jié)合形成彎道螺旋流[25].縱向水流在彎道進口處,主流靠近凸岸側(cè),沿著彎道,主流逐漸向彎道凹岸側(cè)移動,在出口處,主流貼近凹岸.Tang 和Knight[26]將SKM 方法應用于彎道中,基于彎道水流充分發(fā)展的條件下,認為二次流橫向變化率與橫向位置為線性關(guān)系,推導出彎道Ud的解析式.天然彎曲河流常為多個彎道相連,床面形態(tài)存在著深槽和邊灘分布[27],本文嘗試以SKM 方法為基礎(chǔ),推導出彎道中Ud橫向分布計算公式,運用實測資料率定計算參數(shù),將其應用于帶灘槽地形的反向連續(xù)彎道中Ud橫向分布解析計算研究,為下一步的彎道縱向流速快速預測奠定基礎(chǔ).

1 帶灘槽地形彎道垂線平均縱向流速解析解推導

1.1 理論基礎(chǔ)

結(jié)合水流連續(xù)方程,曲線坐標系下,單位體積上的水流縱向動量方程可簡化為

式中,s,r,z分別表示水流的切向、徑向和垂向三個方向,U,V,W和u′,v′,w′分別為s,r,z方向的時間平均流速分量和脈動流速分量,上劃線表示時間平均值,p為壓強,fs為s方向的單位體積力,ρ 為水流密度,g為重力加速度

式中,S0為床面縱比降,S0=??zb/?s,zb為床面高程,H為水深,τrs和τzs分別表示垂直于r和z平面上的切向雷諾應力.

將式(2)～式(4)代入式(1)中,動量方程可表達為

對式(5)沿垂線進行積分可得

式中,τb為床面切應力;1/Sr為床面橫比降;腳標d表示垂線平均值;Se為總水頭能量的比降,用于表示??(zb+H+U2/g)/?s[24];Ua為斷面平均流速;n為糙率,若實驗水槽為水泥面,查表得糙率值為0.014[28];R為水力半徑.

1.2 解析解推導

Tang 和Knight[26]根據(jù)van Balen 等[29]180?彎道試驗數(shù)據(jù),分析了二次流項和二次流附加項的橫斷面分布規(guī)律,發(fā)現(xiàn)其沿橫斷面為線性分布,本文對式(6)右邊二次流項和彎道附加應力項也采用線性分布假設(shè),即

式中,m和k分別為線性分布規(guī)律常數(shù)項和一次項系數(shù),根據(jù)等式左邊可知,m的單位為kg/(m·s2),k的單位為kg/(m2·s2),隨彎道平面形態(tài)及研究斷面所處位置的不同而變化,可采用實測資料率定獲得.

式(6)可寫為

式中,λ 為無量綱渦黏性系數(shù),其取值在河道內(nèi)變化不大,取值為0.07[10,19];f為阻力系數(shù).

將式(13)和式(14)代入式(12)中,可得

本文將帶灘槽地形的橫斷面形態(tài)概化成如圖1所示(圖中灰色區(qū)域為灘地),河槽區(qū)床面高程不變,河灘區(qū)床面線性抬高.對不同的床面形態(tài),由式(16)可求得Ud橫向分布的解析解,解析解的表達式如下.

(a)對河槽區(qū),即水深恒定的區(qū)域,Ud的解析解為

圖1 帶灘槽地形斷面形態(tài)概化圖Fig.1 A sketch of a cross-section with a pool-point bar

(b)對河灘區(qū),即床面線性變化的區(qū)域(床面橫比降為1/Sr),Ud的解析解為

式中,A1～A4為未知的系數(shù),ξ 表示床面上的水深,ξ=H?r/Sr.

1.3 邊界條件

對于Ud的解析式,存在兩種邊界條件,一是無滑移邊界,二是交界面處要滿足連續(xù)性,即:

(1)兩岸(左岸和右岸)邊界處流速為0;

(2)灘槽交界處,主槽和灘地流速相等Ud,i=Ud,i+1(i和i+1 表示參數(shù)取值不同的兩個相鄰區(qū)域);

(3)灘槽交界處,主槽和灘地流速梯度相等,即?Ud,i/?r=?Ud,i+1/?r.

根據(jù)上述邊界條件,解析式中系數(shù)A1,A2,A3,A4均可求解.

2 解析解模型的應用

2.1 帶灘槽地形的連續(xù)彎道縱向流速斷面橫向分布試驗介紹

試驗在武漢大學泥沙重點實驗室的反向連續(xù)彎道水槽中進行,水槽進出口均有一段長7 m 的順直段,上下游兩彎道曲率保持恒定,上游彎道的中心角為90?,下游彎道的中心角為135?,水槽寬度為1.2 m,彎道中心線曲率半徑為2.4 m,兩彎道間的順直連接段長2 m,水槽橫斷面形狀如圖2(b)所示,過渡段中間斷面和進出口5 m 內(nèi)順直段均為矩形斷面,底部高程為0.1 m,彎道橫斷面形狀為凸岸高凹岸低,由凸岸到凹岸,床面高程以1:8 的橫比降由0.2 m 減小為0.1 m(水深減小0.1 m),接著以1:2 的橫比降床面高程減小到0 m (水深減小0.1 m),之后床面高程和水深保持不變,順直段內(nèi)矩形斷面與彎道內(nèi)灘槽地形斷面通過漸變的方式進行連接.

試驗采用三維超聲多普勒流速儀(ADV)進行三維流速測量,儀器采樣頻率為200 Hz,斷面的測點布置橫向間距為5 cm,垂向間距為2 cm,靠近床面區(qū)域加密測量,測點垂向間距減小為1 cm,由于測量探頭的大小及儀器采樣體積為探頭下方5 cm,因此在邊壁和水面5 cm 以內(nèi)無法進行測量,縱向流速的垂線平均值根據(jù)各測線的流速數(shù)據(jù)取垂線平均獲得.

圖2 連續(xù)彎道試驗布置Fig.2 Layout of the consecutive bends experiment

三組工況的水力學條件如表1 所示,典型天然彎曲河流的Fr范圍為0.07～0.26[30],三種試驗工況的Fr取值在此范圍內(nèi),并保持流量不變,研究水深變化的影響,采用試驗測量的五個斷面(斷面A～E)進行驗證.表中,Q為進口流量,H為出口斷面水深,U為進口斷面平均流速,Fr為佛汝德數(shù),Re為雷諾數(shù).

表1 試驗條件Table 1 Experimental conditions

2.2 結(jié)果與分析

圖3 為解析解的擬合曲線與試驗結(jié)果的比較(圖中灰色區(qū)域為灘地),圖中圓點為根據(jù)實測的三維流速數(shù)據(jù)計算出的Ud,由圖可知實測的Ud沿橫斷面的分布規(guī)律為深槽流速較小,流速沿橫向變化率較小,灘地上流速較大,流速沿橫向變化率增大,Ud的最大值在靠近凸岸側(cè),下游彎道中Ud橫向分布較上游彎道更不均勻,由工況一到工況三,出口水深逐漸減小即水流流速增大,在下游彎道中Ud的最大值逐漸向凹岸偏移.

圖3 流速計算值與實測值對比Fig.3 Comparison of modelled depth-averaged velocity distribution with experimental data

圖3 中曲線為計算出的擬合曲線,水深恒定不變的河槽區(qū)(距凹岸0～0.2 m 范圍內(nèi))采用式(17)進行擬合,床面線性變化的邊坡區(qū)(距凹岸0.2～1.2 m 范圍內(nèi))由式(20)擬合,參數(shù)m的取值根據(jù)地形分為三個區(qū)域,分別為水平段m1、橫比降為1/2 的斜坡段m2和橫比降為1/8 的斜坡段m3,參數(shù)k的取值分為兩個區(qū)域(靠近凹岸側(cè)為k1,靠近凸岸側(cè)為k2),分區(qū)位置在主流附近,解析解表達式中的系數(shù)A1～A4可根據(jù)邊界條件求解.由圖可知擬合曲線與實測數(shù)據(jù)吻合良好,在床面地形變化處,Ud沿程變化率較大時以及Ud峰值的大小與位置,計算值對實測值均能較好擬合,隨著出口水深減小,Ud最大值的位置變化,計算結(jié)果也能較好地反映這一規(guī)律,僅在靠近邊壁處擬合曲線與實測數(shù)據(jù)稍有偏差,凹岸側(cè)偏差可能是凹岸次生反向小環(huán)流影響導致,凸岸邊壁附近偏差則可能是由于水深較淺,水流結(jié)構(gòu)較為復雜導致,各斷面的均方根誤差均小于0.015 m/s,平均相對誤差均在5%以內(nèi),實測值與計算值之間誤差較小.綜上所述,在帶有灘槽的反向連續(xù)彎道中,此解析式能較好擬合不同水流條件下沿程各斷面垂線平均縱向流速的橫向分布規(guī)律.

3 計算參數(shù)敏感性分析及取值

對參數(shù)m和k的取值進行敏感性分析,此處以工況二的斷面A擬合結(jié)果為基礎(chǔ)進行分析.圖4對k的數(shù)值大小進行了敏感性分析,在0～1.1 m 和1.1～1.2 m 取值不同(分別表示為k1和k2),將k1和k2的取值分別在原基礎(chǔ)上增大和減小0.05 (每次僅變化一個參數(shù)的取值),由圖可知k1數(shù)值變化對凹岸深槽區(qū)域流速計算值基本無影響,主要變化發(fā)生在靠近凸岸的邊灘附近,且計算值變化較小,k2的變化僅影響1～1.2 m 內(nèi)的局部區(qū)域,流速計算值變化較小,流速最大值的大小和位置未發(fā)生明顯變化;圖5 對k1和k2分區(qū)位置的敏感性進行了分析,虛線Lk表示分界位置,將初始分界位置分別向左和向右移動0.1 m,由圖可以發(fā)現(xiàn),分界位置的變化對流速的峰值影響較大,流速最大值變化明顯,且位置也發(fā)生變化,由此可將主流線的位置作為k分區(qū)的初步依據(jù).

圖4 k 變化時垂線平均流速計算值Fig.4 Modelled depth-averaged velocity for various k values

圖5 k 分區(qū)調(diào)整時垂線平均流速計算值Fig.5 Modelled depth-averaged velocity for various k values in different panels

參數(shù)m分為三個區(qū)域進行取值,為m1(0～0.2 m),m2(0.2～0.4 m)和m3(0.4～1.2 m),圖6 為對參數(shù)m取值的敏感性分析結(jié)果,分別對m1,m2和m3的初始值增大和減小0.05 (每次僅變化一個參數(shù)的取值),由圖可知,m1的取值增大或減小對流速計算值基本無影響,m2的取值變化僅在0.2～0.4 m 附近的流速計算值產(chǎn)生影響,流速計算值變化較小,m3取值變化對整個斷面的流速計算值基本均有影響,在0.4～1.2 m 區(qū)域內(nèi),流速計算值發(fā)生明顯變化;圖7為m分區(qū)位置變化的敏感性分析結(jié)果,圖中虛線Lma和Lmb分別為m1與m2和m2與m3分區(qū)的位置,將初始分界位置分別向左和向右移動0.1 m(每次僅變化一個分界位置),由圖可知m1和m2分界位置向斜坡段移動對流速的計算值影響較為明顯,m2和m3分界位置的變化對流速計算值影響較小.

圖6 m 變化時垂線平均流速計算值Fig.6 Modelled depth-averaged velocity for various m values in different panels

圖7 m 分區(qū)位置變化時垂線平均流速計算值Fig.7 Modelled depth-averaged velocity for various panels of m

參數(shù)m和k的取值如表2 所示,Lk1和Lk2分別為k1和k2取值區(qū)域的長度,由表中數(shù)據(jù)可知,m和k的值與Tang 和Knight[26]文中的取值是相近的.沿斷面水平段m1的取值小于斜坡段m2和m3的取值,可能與凹岸主槽側(cè)流速較小有關(guān),參數(shù)m的取值整體上隨水深的減小逐漸減小,水深減小到一定程度,河灘區(qū)m2和m3的取值由正變?yōu)樨?k1和k2取值分區(qū)的界限與主流位置較為接近,分區(qū)的位置變化較小,可能與主流位置變化較小相關(guān),k1和k2取值的相對大小整體上由工況一的k1大于k2變?yōu)楣r三的k1小于k2.根據(jù)參數(shù)的取值敏感性分析可知m1,m2和k1,k2的取值變化對計算結(jié)果的影響較小,計算結(jié)果對m3的值較為敏感,但m3的取值沿程變化較小,因此給出各參數(shù)沿水槽的均值作為初步取值參考.

4 討論

本文對二次流項和二次流附加項采用與Tang 和Knight[26]相同的線性假設(shè),其僅分析彎道出口附近135?斷面的數(shù)據(jù),本文將解析解用于彎道中的多個斷面,此處增加繪制van Balen 等[29]的180?彎道進口附近29?斷面的二次流項ρH(UV)d和二次流附加項的斷面分布(如圖8 和圖9所示),由圖可知其在彎道進口斷面也符合線性假設(shè).進一步分析二次流項和二次流附加項在本文試驗水槽中的斷面分布規(guī)律,選取工況二csC斷面實測數(shù)據(jù)進行討論,由圖10 和圖11 可知二次流項和二次流附加項的斷面分布符合線性分布規(guī)律,由此可知線性假設(shè)在本試驗中也是符合的.由此可知,對二次流項和二次流附加項沿橫斷面的線性分布假設(shè)在彎道沿程多個斷面均適用,并在本文連續(xù)彎道中也是符合的,表明本文采取的線性假設(shè)在本文的研究問題中是合理的.

表2 參數(shù)取值Table 2 The values of parameters

圖8 二次流項ρH(UV)d29?橫斷面分布[29]Fig.8 Lateral distribution of secondary flow term,ρH(UV)dat cross-section of 29?[29]

圖9 二次流附加項dz 29?橫斷分布[29]Fig.9 Lateral distribution of additional secondary flow term,dz at cross-section of 29?[29]

圖10 二次流項ρH(UV)d橫斷面分布(csC)Fig.10 Lateral distribution of secondary flow term,ρH(UV)d(csC)

本文推導出彎道中的Ud沿橫斷面解析解表達式,與帶灘槽地形的反向連續(xù)彎道實測數(shù)據(jù)符合良好,根據(jù)實測數(shù)據(jù)率定的參數(shù)m和k的取值與Tang和Knight[26]的取值是相近的,表明參數(shù)的取值對與不同的彎道存在相近的范圍,并分析了參數(shù)取值隨水位變化的規(guī)律,但具體的參數(shù)取值與水流條件和彎道形態(tài)之間的關(guān)系仍需更加深入的研究,以更好地推廣應用解析解公式.

圖11 二次流附加項橫斷分布(csC)Fig.11 Lateral distribution of additional secondary flow term,

5 結(jié)論

本文基于SKM 方法,對動量方程進行垂線積分,推導出了彎道中Ud的橫向分布解析解公式,將其應用于帶灘槽地形的反向連續(xù)彎道中.

(1)將解析解公式應用于有灘槽的反向連續(xù)彎道,調(diào)整線性假設(shè)中參數(shù)的值,可以使計算值與實測值在不同的出口水流條件下,對沿程各斷面流速值橫向分布規(guī)律都能吻合良好,表明本文的線性假設(shè)具有一定的合理性.

(2)參數(shù)取值,一次項系數(shù)k的大小對流速橫向分布趨勢有較大影響,k的分區(qū)對主流有較大的影響,調(diào)整分區(qū)可使計算值更好擬合灘地區(qū)流速的峰值大小和位置.常數(shù)項m整體上隨水深的減小而減小,可對m按照斜率變化分區(qū)對局部流速計算值進行調(diào)整,河槽區(qū)m的取值小于河灘區(qū),流速計算值對m在河槽區(qū)與河灘區(qū)間分界線的變化較為敏感,河灘區(qū)內(nèi)m的分區(qū)則對流速計算值影響較小.

(3)本文推導的解析解公式應用范圍較廣,可用于地形變化的反向連續(xù)彎道中,但對于該模式的應用需要進行更加詳細的研究:①增加試驗工況,進一步分析參數(shù)m和k取值與進出口水流條件(流量、水位和斷面平均流速等)和彎道形態(tài)(曲率、比降和過渡段長度等)之間的規(guī)律,總結(jié)參數(shù)分區(qū)依據(jù)和合理的取值范圍;②分析m和k的物理意義,利用更加詳細精確的數(shù)據(jù)計算其表示的各項具體分布規(guī)律,以更深入理解m和k取值變化規(guī)律的具體原因.