李佳佳，李志偉，張長寬，陳永平，3

(1.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇南京 210098;2.香港理工大學(xué)土木及環(huán)境工程學(xué)系，香港; 3.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇南京 210098)

在明渠湍流的研究中，人們通過具有光滑或粗糙床面的湍流模型來認(rèn)識其流速、雷諾應(yīng)力和紊動強(qiáng)度的垂線分布等復(fù)雜水力特性，這些水力特性對于明渠湍流的床面切應(yīng)力、水流中泥沙起動和輸移以及污染物質(zhì)擴(kuò)散等的研究都具有十分重要的意義。針對光滑床面明渠水流的模擬技術(shù)已較為成熟，但光滑床面是一種理想化的狀態(tài)，在實際河道中，床面粗糙不平，復(fù)雜多變，因此研究粗糙床面明渠湍流水力特性更有助于認(rèn)識天然河道的實際情況。

常見的計算流體力學(xué)模型多基于RANS方程(雷諾平均NS方程)，然而在實際水流中，尤其是近床區(qū)域，存在許多小尺度渦旋，使得時均流具有強(qiáng)烈的空間不均勻性，因而RANS方程并不適用。為了解決這一問題，有學(xué)者［1-4］提出一種基于DANS方程(雙向平均NS方程)的模擬方法，先后在時間和空間上對NS方程進(jìn)行平均，即在傳統(tǒng)RANS方程的基礎(chǔ)上再進(jìn)行空間平均。該方法最初應(yīng)用于氣象研究［5］，后來逐漸被運用到水力學(xué)研究領(lǐng)域［1-4］。經(jīng)過空間平均之后，DANS方程出現(xiàn)了一些新的物理項，如形狀應(yīng)力和形狀阻力等。一般來講，形狀應(yīng)力比雷諾應(yīng)力小得多，在以往的研究中通常被忽略，但是近期的研究表明:形狀應(yīng)力并不是在任何情況下都可以被忽略的。Poggi等［6-7］用物理試驗來研究植物分布密度對水力特性的影響，結(jié)果表明，當(dāng)?shù)撞恐参锓植嫉拿芏群芨邥r，形狀應(yīng)力可以被忽略;當(dāng)?shù)撞恐参锓植嫉帽容^稀疏時，形狀應(yīng)力在植被以上的區(qū)域可以被忽略，而在植被分布層內(nèi)的區(qū)域不能被忽略。Mclean等［8］用物理試驗來模擬河床底部為沙丘時的明渠水力特性，給出了沙丘以上水域的試驗數(shù)據(jù)，結(jié)果表明，在沙丘以上水域的形狀應(yīng)力可以忽略不計。Campbell等［9-10］用物理試驗研究了礫石河床的相關(guān)狀況，研究表明:形狀應(yīng)力在礫石頂部以上區(qū)域可以忽略，而在其以下區(qū)域則不能忽略，其貢獻(xiàn)程度與雷諾應(yīng)力相當(dāng)。

從以上的研究成果可以看出，河床底部障礙物的材質(zhì)及分布的疏密情況均會給明渠水力特性帶來影響。Nikora等［3］研究發(fā)現(xiàn)，障礙物垂向分布的不均勻性也同樣會給明渠的水力特性帶來一定的影響，他將障礙物的垂向分布情況用孔隙率φ來表示，φ被定義為水的體積與總體積的比值，用以反映河床底部的地形分布狀況。Nikora等［1-2］還系統(tǒng)地給出了考慮φ之后DANS方程推導(dǎo)的理論基礎(chǔ)，并將其應(yīng)用于明渠水流的模擬研究中。到目前為止，φ的相關(guān)研究成果并不多。

本文在前人研究的基礎(chǔ)上，基于DANS方程建立粗糙面明渠水流數(shù)學(xué)模型，研究粗糙床面明渠水流的形狀應(yīng)力，定量探討其存在對雷諾應(yīng)力和速度分布的影響，并研究障礙物孔隙率的分布對形狀應(yīng)力、雷諾應(yīng)力及垂線速度分布的影響。

1 數(shù)值模型的建立

為了簡化模型，本文僅考慮均勻流的情形，DANS方程可以簡化為一維的形式，忽略二次流影響后，其連續(xù)性方程和動量方程表述如下。

連續(xù)性方程:

動量方程:

S-A湍流模型是從經(jīng)驗和量綱分析出發(fā)［11］，由針對簡單流動逐漸發(fā)展為適用于帶有層流流動的固壁湍流流動的方程模型。

相對于兩方程湍流模型，S-A湍流模型的計算量較小，穩(wěn)定性好，計算網(wǎng)格在物面處不需要特別精細(xì)。在S-A湍流模型中，雷諾應(yīng)力可表示為

于是式(2)變成:

式中:υt——湍流黏性系數(shù)，需要用S-A模型進(jìn)行計算。

S-A湍流模型控制方程如下:

式中:L——流場中某點至壁面的最小距離;υ——S-A湍流模型中的黏性系數(shù)，通過S-A湍流模型控制方程給出;wz——旋轉(zhuǎn)角速度。

與傳統(tǒng)的RANS方程相比，DANS方程(式(2))多了新的一項該項的處理可以借鑒LES(大渦數(shù)值模擬)的思想［12-14］，即將形狀應(yīng)力視為LES內(nèi)的亞格子應(yīng)力項，采用與LES類似的亞格子渦黏模型來計算形狀應(yīng)力項:

式中:υs——亞格子渦黏系數(shù);Cs——經(jīng)驗系數(shù);Δ——網(wǎng)格空間步長。

將式(7)代入式(4)，DANS控制方程可改寫為

與傳統(tǒng)的RANS方程相比，DANS方程還考慮了障礙物φ的影響［15］。對于障礙物以上的純水區(qū)域，φ= 1;對于障礙物以下的區(qū)域，0≤φmin≤φ＜1。本文針對不同孔隙率分布狀況對水流結(jié)構(gòu)的影響進(jìn)行研究對比。

2 形狀應(yīng)力對雷諾應(yīng)力和垂向速度分布的影響

首先定量探討形狀應(yīng)力的大小、分布及其對雷諾應(yīng)力和速度分布的影響。計算值來源于本文建立的基于DANS方程的粗糙床面明渠水流數(shù)學(xué)模型模擬結(jié)果，用于驗證的實測資料來源于Nikora等［3］的物理試驗。該試驗在一個長為12 m、寬為0.75 m的水槽內(nèi)進(jìn)行的，水槽底部被直徑為63.9 mm的球體覆蓋，采用聲學(xué)多普勒流速儀(ADV)進(jìn)行測量，分別測量了2種工況下的數(shù)據(jù)，見表1。表1中，Q為流量;H為水深;d0為障礙物高度;Re為雷諾數(shù);u*為摩阻流速，反映壁面處阻力的大小，用于后續(xù)數(shù)據(jù)的無量綱化過程。

表1 文獻(xiàn)［3］試驗參數(shù)設(shè)置Table 1 Key parameters in experiment of literature［3］

圖1分別為2種工況下粗糙床面明渠水流數(shù)學(xué)模型在Cs為0.05、0.10、0.15時形狀應(yīng)力、雷諾應(yīng)力以及速度分布計算值，并與文獻(xiàn)［3］的物理試驗測量值進(jìn)行對比驗證。圖中縱坐標(biāo)由下至上代表從床底到水面的過程，0～1為障礙物分布層內(nèi)，坐標(biāo)1以上代表障礙物頂部以上至水面的水流區(qū)域。

從圖1(a)(b)可看出，形狀應(yīng)力在障礙物分布層內(nèi)是不可忽略的，而在障礙物上部則逐漸減小直至可忽略不計。由圖1(c)(d)可見，Cs=0時，形狀應(yīng)力值為零，雷諾應(yīng)力值最大。當(dāng)z/d0=1時，即在障礙物頂部的位置時，雷諾應(yīng)力達(dá)到最大值，并向水面線性減小，向底部也是逐漸減小的，計算結(jié)果與文獻(xiàn)［3］的試驗結(jié)果一致。雷諾應(yīng)力由水面至障礙物頂部的線性增加是由重力分量的作用引起，而其從障礙物頂部至河床逐漸減小則是由于在障礙物層內(nèi)存在拖曳力及形狀應(yīng)力與其共同保持與重力分量的平衡。圖1(e)(f)為速度的垂向分布圖，由計算值與實測值對比可知，在障礙物分布層以內(nèi)，速度分布呈線性分布，而在其上部則符合對數(shù)分布規(guī)律。速度值隨形狀應(yīng)力的增大而減小，這與理論預(yù)期是一致的，當(dāng)水流中引入的阻力越大，水流速度就相應(yīng)變得越小。

圖1 不同Cs條件下形狀應(yīng)力、雷諾應(yīng)力以及速度的垂向分布Fig.1 Vertical distributions form-induced stress，Reynolds stress，and velocity with different values of Cs

3 障礙物孔隙率對水流特性的影響

同樣采用文獻(xiàn)［3］的試驗進(jìn)行障礙物孔隙率分布的影響研究，孔隙率分布狀況見圖2。0.5常數(shù)分布是指在障礙物層內(nèi)φ為常數(shù)0.5，在障礙物層以上即純水區(qū)域為1;0.5線性分布是指在河床底部障礙物φ為0.5，并且在障礙物層從河床底部到障礙物頂部由0.5線性變化至1;同樣，0.5分段分布指的是在河床底部至障礙物層一半高度處φ為0.5，并從障礙物一半高度處由0.5至1線性變化至障礙物頂部純水區(qū)域;此外在研究過程中還取了1.0常數(shù)分布，指的是在計算域內(nèi)忽略障礙物分布的影響。計算值與實測值的對比驗證結(jié)果見圖3。

圖3為在不同的孔隙率分布狀況下形狀應(yīng)力、雷諾應(yīng)力以及速度的垂向分布圖。從圖3(a)(b)可以看出，不同的φ分布狀況對形狀應(yīng)力分布有一定的影響，適當(dāng)改變φ分布，計算值會更加接近實測值。從圖3(c)(d)可以看出，當(dāng)φ為0.5常數(shù)分布時，由于φ在障礙物頂端處存在突變，因而雷諾應(yīng)力在此處也有一定的突變。從圖3(e)(f)中看出，φ的分布對其速度分布影響并不大，并且φ分布越光滑對其影響越小。

圖2 孔隙率φ的分布Fig.2 Distribution of porosity φ

圖3 不同孔隙率分布條件下形狀應(yīng)力、雷諾應(yīng)力以及速度的垂向分布Fig.3 Vertical distributions of form-induced stress，Reynolds stress，and velocity with different distributions of porosity

4 結(jié) 語

本文基于DANS方程，運用修正的S-A湍流模型，采用渦黏系數(shù)法建立了粗糙床面明渠水流數(shù)學(xué)模型，用該模型計算了速度、雷諾應(yīng)力及形狀應(yīng)力的大小及分布并與實測值進(jìn)行比較，此外還研究了底部障礙物的孔隙率分布情況對形狀應(yīng)力、雷諾應(yīng)力以及垂線速度分布的影響，得到以下主要結(jié)論:

a.形狀應(yīng)力在障礙物分布層以內(nèi)是不可忽略的，而在障礙物上部會迅速減小直至忽略不計。

b.形狀應(yīng)力的存在造成障礙物頂部附近的雷諾應(yīng)力明顯減小，但在障礙物頂部以上至水面區(qū)域，形狀應(yīng)力對雷諾應(yīng)力的分布幾乎無影響。

c.形狀應(yīng)力的存在使得速度有一定程度的減小，并且形狀應(yīng)力值越大，速度會隨之變得越小。

d.障礙物孔隙率的分布對速度、雷諾應(yīng)力及形狀應(yīng)力等會存在一定的影響，并且當(dāng)孔隙率的分布越光滑時，其影響越小。

在實際情況中，由于地形分布十分復(fù)雜，孔隙率的分布只能用理想化模型代替，由于目前有關(guān)孔隙率分布的實測驗證資料較少，在數(shù)學(xué)模型研究中采用哪一種分布狀況更為合理還有待進(jìn)一步研究。

［1］NIKORA V，MCEWAN I，MCLEAN S，et al.Double-averaging concept for rough-bed open-channel and overland flows: theoretical background［J］.ASCE，Journal of Hydraulic Engineering，2007，133(8):873-883.

［2］NIKORA V，MCLEAN S，COLEMAN S，et al.Double-averaging concept for rough-bed open-channel and overland flows: applications［J］.ASCE，Journal of Hydraulic Engineering，2007，133(8):884-895.

［3］NIKORA V，GORING D，MCEWAN I，et al.Spatially averaged open-channel flow over rough bed［J］.ASCE，Journal of Hydraulic Engineering，2001，127(2):123-133.

［4］SMITH J D，MCLEA S R.Spatially averaged flow over a wavy surface［J］.Journal of Geophysical Research，1977，82(12): 1735-1746.

［5］WILSON N R，SHAW R H.A higher order closure model for canopy flow［J］.Journal of Applied Meteorology and Climatology， 1977，16:1197-1205.

［6］POGGI D，KATUL G G，ALBERTSON J D.A note on the contribution of dispersive fluxes to momentum transfer within canopies［J］.Boundary-Layer Meteorology，2004，111:615-621.

［7］POGGI D，PORPORATO A，RIDOLFI L，et al.The effect of vegetation density on canopy sub-layer turbulence［J］.Boundary-Layer Meteorology，2004，111:565-587.

［8］MCLEAN S，NIKORA V.Characteristics of turbulent unidirectional flow over rough beds:double-averaging perspective with particular focus on sand dunes and gravel beds［J］.Water Resources Research，2006，42，W10409，DOI:10.1029/ 2005WR004708.

［9］CAMPBELL L，MCEWAN I，NIKORA V，et al.Bed load effects on hydrodynamics of rough bed open channel flows［J］.ASCE，Journal of Hydraulic Engineering，2005，131(7):576-585.

［10］ABERLE J.Spatially averaged near-bed flow field over rough armor layers［C］//ELSA C T，CARDOSO A H，LEAL J G，et al. Proceedings of the International Conference on Fluvial Hydraulics.Lisbon，Portugal:Taylor＆Francis，2006:153-162.

［11］SPALART P，ALLMARAS R.A one-equation turbulence model for aerodynamic flows［C］//IANNOTTA B.30th Aerospace Sciences Meeting and Exhibit.Reno，USA:AIAA，1992:439.

［12］ZENG C，LI C W.Modeling flows over gravel beds by a drag force method and a modified S-A turbulence closure［J］.Advances in Water Resources，2012，46:84-95.

［13］LI C W，YAN K.Numerical investigation of wave-current-vegetation interaction［J］.ASCE，Journal of Hydraulic Engineering，2007，133(7):794-803.

［14］LI C W，XIE J F.Numerical modeling of free surface flow over submerged and highly flexible vegetation［J］.Advances in Water Resources，2011，34:468-477.

［15］MANES C，POKRAJAC D，MCEWAN I，et al.Turbulence structure of open channel flows over permeable and impermeable beds: a comparative study［J］.Physics of Fluids，2009，21(12):5109.