基于DANS方程的粗糙床面明渠水力特性研究

李佳佳，李志偉，張長(zhǎng)寬，陳永平，3

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

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

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

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

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

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

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

連續(xù)性方程:

動(dòng)量方程:

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

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

于是式(2)變成:

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

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

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

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

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

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

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

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

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

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

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

從圖1(a)(b)可看出，形狀應(yīng)力在障礙物分布層內(nèi)是不可忽略的，而在障礙物上部則逐漸減小直至可忽略不計(jì)。由圖1(c)(d)可見，Cs=0時(shí)，形狀應(yīng)力值為零，雷諾應(yīng)力值最大。當(dāng)z/d0=1時(shí)，即在障礙物頂部的位置時(shí)，雷諾應(yīng)力達(dá)到最大值，并向水面線性減小，向底部也是逐漸減小的，計(jì)算結(jié)果與文獻(xiàn)［3］的試驗(yàn)結(jié)果一致。雷諾應(yīng)力由水面至障礙物頂部的線性增加是由重力分量的作用引起，而其從障礙物頂部至河床逐漸減小則是由于在障礙物層內(nèi)存在拖曳力及形狀應(yīng)力與其共同保持與重力分量的平衡。圖1(e)(f)為速度的垂向分布圖，由計(jì)算值與實(shí)測(cè)值對(duì)比可知，在障礙物分布層以內(nèi)，速度分布呈線性分布，而在其上部則符合對(duì)數(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 障礙物孔隙率對(duì)水流特性的影響

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

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

圖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é) 語(yǔ)

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

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

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

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

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

在實(shí)際情況中，由于地形分布十分復(fù)雜，孔隙率的分布只能用理想化模型代替，由于目前有關(guān)孔隙率分布的實(shí)測(cè)驗(yà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.