白 靜方紅衛(wèi)何國(guó)建謝崇寶高 虹

?(中國(guó)灌溉排水發(fā)展中心，北京100054)

?(清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100084)

??(北京中灌綠源國(guó)際咨詢有限公司，北京100054)

細(xì)顆粒泥沙凈沖刷和輸移的大渦模擬研究

白 靜?方紅衛(wèi)?,1)何國(guó)建?謝崇寶?高 虹??

?(中國(guó)灌溉排水發(fā)展中心，北京100054)

?(清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100084)

??(北京中灌綠源國(guó)際咨詢有限公司，北京100054)

在傳統(tǒng)水沙輸移數(shù)值模擬研究中一般采用雷諾時(shí)均模擬技術(shù)(Reynolds-averaged simulation,RANS).與RANS相比，大渦模擬技術(shù)(large eddy simulation,LES)能夠更加精確反映細(xì)部流動(dòng)結(jié)構(gòu)，計(jì)算機(jī)的發(fā)展使得采用LES探討水流和泥沙運(yùn)動(dòng)規(guī)律成為可能.本文嘗試給出凈沖刷條件下懸沙計(jì)算的邊界條件，采用動(dòng)態(tài)亞格子模式對(duì)循環(huán)槽道和長(zhǎng)槽道中的水流運(yùn)動(dòng)和泥沙輸移進(jìn)行了三維大渦模擬研究.利用直接數(shù)值模擬(direct numerical simulation,DNS)結(jié)果對(duì)LES模型進(jìn)行了率定，計(jì)算結(jié)果符合良好，在此基礎(chǔ)上初步探討了泥沙濃度、湍動(dòng)強(qiáng)度和湍動(dòng)通量等的分布特征.結(jié)果表明，凈沖刷條件下輸沙平衡時(shí)泥沙濃度符合Rouse公式分布，單向流動(dòng)中泥沙濃度沿著流向逐漸增大.泥沙濃度湍動(dòng)強(qiáng)度和湍動(dòng)通量都在近底部達(dá)到最大值，沿著垂向迅速減小.湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)基本上在水深中間處達(dá)到最大.湍動(dòng)Schmidt數(shù)沿著水深方向不是常數(shù)，在近底部和自由水面附近較大，在水深中間處較小.

細(xì)顆粒泥沙,凈沖刷,大渦模擬,泥沙輸移

引言

泥沙輸移一直是水力學(xué)及河流動(dòng)力學(xué)領(lǐng)域研究的熱門問題.淤積和沖刷是不平衡輸沙的重要形式.不平衡輸沙是河床變形和水庫(kù)淤積的根源.細(xì)顆粒泥沙的輸移一直受到廣泛的關(guān)注[1-2].細(xì)顆粒泥沙在水體污染、水體富營(yíng)養(yǎng)化和水生生態(tài)系統(tǒng)健康等方面有重要的影響.來自農(nóng)業(yè)、工業(yè)和生活等污水進(jìn)入河道后，以細(xì)顆粒泥沙為載體的生物膜生長(zhǎng)過程，改變了泥沙的形態(tài)特征、沉降和起動(dòng)特性[3].細(xì)顆粒泥沙的吸附解析作用和絮凝過程，對(duì)河道和水庫(kù)中污染物的遷移轉(zhuǎn)化[4]和重金屬的遷移[5]有復(fù)雜的影響.泥沙輸移的計(jì)算方法主要有兩種：歐拉方法[6-7]和拉格朗日方法，歐拉方法的計(jì)算量大，只能模擬有限泥沙顆粒的運(yùn)動(dòng).在泥沙輸移和河床變形模擬中，泥沙顆粒的數(shù)量巨大，因此應(yīng)當(dāng)采用拉格朗日方法進(jìn)行模擬.目前絕大多數(shù)的泥沙運(yùn)動(dòng)和水流模擬成果主要采用雷諾時(shí)均模擬(Reynolds-averaged simulation,RANS)完成[8-10]，Keylock等[11]認(rèn)為隨著計(jì)算機(jī)技術(shù)水平的提高，使采用大渦模擬(large eddy simulation,LES)追蹤和研究湍流相干結(jié)構(gòu)的變化成為可能，由于湍流相干結(jié)構(gòu)對(duì)泥沙起動(dòng)和湍流輸運(yùn)有極其重要的作用，LES在泥沙輸移和河床演變研究中的潛力和優(yōu)勢(shì)也逐漸表現(xiàn)出來.根據(jù)前人的研究，對(duì)于同樣的算例,LES的計(jì)算網(wǎng)格量為三維雷諾時(shí)均模擬的數(shù)十倍乃至數(shù)百倍[12]，計(jì)算時(shí)間和計(jì)算量隨著計(jì)算區(qū)域長(zhǎng)度的增大而增加.

國(guó)內(nèi)外的學(xué)者針對(duì)LES中泥沙計(jì)算邊界、泥沙起動(dòng)的機(jī)理、濃度分布特征和沙紋的形成過程進(jìn)行了探討.在定床計(jì)算中，泥沙顆粒一般較細(xì)，一般只涉及到懸移質(zhì).Zedler等[13-14]假定底部泥沙濃度的梯度等于pick-up函數(shù)與擴(kuò)散系數(shù) ε的商，即?s/?z=-Pk/ε，先后模擬了單向流和振蕩流中泥沙起動(dòng)過程，在單向流中得到了泥沙濃度隨時(shí)間的變化過程，得到了振蕩流中流速、壓強(qiáng)以及泥沙濃度分布.其研究認(rèn)為泥沙的起動(dòng)原因有兩種：一種是泥沙在渦體形成過程中被帶起，一種由于床面切應(yīng)力較高，泥沙直接起動(dòng).Chou等[15]利用動(dòng)態(tài)混合模型和LES計(jì)算了泥沙在明渠中起動(dòng)直至平衡的全過程，平衡時(shí)的泥沙濃度與Rouse公式符合較好，認(rèn)為通過床面進(jìn)入水體的泥沙通量主要由流向渦決定.Widera等[16-17]利用LES計(jì)算平衡時(shí)泥沙在平底和有沙波的條件下槽道水流中的泥沙分布，主要探討了湍動(dòng)Schimidt數(shù)的取值問題，指出了在RANS中將不同水深處湍動(dòng)Schmidt數(shù)取為常數(shù)的做法是不可靠的；在有沙波存在的明渠中，底部泥沙濃度受地形的影響較大，在遠(yuǎn)離床底處，泥沙濃度與Rouse公式符合很好.靠近床面處，湍動(dòng)Schmidt數(shù)對(duì)地形不是很敏感.

在前面所提及的利用LES對(duì)泥沙分布與輸移的研究中，Zelder等[13-14]采用的泥沙邊界條件，在擴(kuò)散系數(shù)ε較小時(shí)，容易造成泥沙濃度s在底部的梯度過大.Chou等[15]引入泥沙沉降通量Db和起動(dòng)通量Eb的概念，Eb通過pick-up函數(shù)計(jì)算得到，近底處的泥沙濃度sb通過豎直方向上鄰近的三點(diǎn)濃度插值得到，這種泥沙邊界的處理方法，在計(jì)算時(shí)沒有用到參考高度的概念，比較簡(jiǎn)單，但是Eb實(shí)際上是全沙通量，并不準(zhǔn)確.由于流向采用周期性邊界條件，只能給出泥沙濃度和床面變形隨時(shí)間的變化過程，不能給出空間變化過程.Widera等[16-17]在水平方向上采用周期性邊界條件，在底部泥沙計(jì)算中采用零通量的邊界條件，使得計(jì)算結(jié)果只能反映泥沙輸沙平衡后的最終狀態(tài)，不能反映泥沙隨著時(shí)間和空間的變化情況.

泥沙可以分為懸移質(zhì)和推移質(zhì)，推移質(zhì)集中在床面附近.與泥沙的凈淤積邊界條件相比，凈沖刷條件相對(duì)比較復(fù)雜[18].為了解決凈沖刷過程中懸沙的起算高度問題，在計(jì)算中保證泥沙和水流計(jì)算網(wǎng)格的一致性，本文中將借鑒RANS的處理方式，在懸沙計(jì)算中引入?yún)⒖几叨雀拍睿⒓?xì)顆粒泥沙凈沖刷條件下的近底邊界條件，嘗試給出輸沙平衡條件下和單向流動(dòng)中的水流運(yùn)動(dòng)和泥沙輸移規(guī)律，對(duì)其中水流和泥沙的時(shí)均特征、湍動(dòng)特性進(jìn)行分析和討論，進(jìn)一步統(tǒng)計(jì)湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)的分布，給出湍動(dòng)Schmidt數(shù)的分布規(guī)律.

1 控制方程

水流計(jì)算的控制方程為

其中“ˉ”表示過濾后的變量，ui(i=1,2,3)表示i方向上的無量綱速度；υ是雷諾數(shù)Re的倒數(shù)；p是無量綱壓強(qiáng)；對(duì)控制方程中的非線性對(duì)流項(xiàng)過濾產(chǎn)生的亞格子應(yīng)力為它代表亞格子尺度運(yùn)動(dòng)對(duì)可解尺度運(yùn)動(dòng)的影響.

其中υSGS為亞格子黏性系數(shù)為 Kronecker算子，本文采用Germano等[19]提出的動(dòng)態(tài)模式進(jìn)行求解亞格子黏性系數(shù)υSGS.

懸沙運(yùn)動(dòng)的控制方程為

其中s為無量綱的懸沙濃度，ε為懸沙的分子擴(kuò)散系數(shù)，ws為懸沙的無量綱沉速，ε=υ/σ，σ是Schmidt數(shù)，本文中取代表了亞格子尺度的水流運(yùn)動(dòng)對(duì)懸沙的輸運(yùn)作用，?i為L(zhǎng)ES中懸沙控制方程的不封閉項(xiàng).

其中εSGS為泥沙亞格子尺度的擴(kuò)散系數(shù)，εSGS與亞格子尺度黏性系數(shù)有關(guān)：εSGS=υSGS/σSGS，σSGS是亞格子尺度的Schmidt數(shù)，本文中取1.

2 數(shù)值方法與懸沙參考高度

對(duì)德國(guó)卡爾斯魯厄大學(xué)流體力學(xué)研究所開發(fā)LESOCC2程序[20-23]進(jìn)行了修改，加入了懸沙計(jì)算模塊.在程序中采用有限體積法對(duì)方程進(jìn)行離散，變量布置在非交錯(cuò)曲線網(wǎng)格上，分別采用二階精度的中央差分和HLPA差分[24]近似水流和泥沙計(jì)算控制方程中的對(duì)流通量，統(tǒng)一采用中央差分近似擴(kuò)散通量，時(shí)間推進(jìn)利用二階精度的3步Runge--Kutta法.在Runge--Kutta法的第3步，采用Stone[25]提出的強(qiáng)隱式迭代(strongly implicit procedure,SIP)方法求解壓力修正的泊松方程.

推移質(zhì)集中在床面表面附近區(qū)域，懸移質(zhì)計(jì)算的網(wǎng)格起點(diǎn)也就是參考高度a，應(yīng)該不低于推移質(zhì)層的高度.不同的研究學(xué)者對(duì)參考高度有著不同的研究，Celik等[26]認(rèn)為：對(duì)于沒有沙波的槽道，參考高度為a=2/3ks=2d50，其中ks為床面的當(dāng)量粗糙高度；對(duì)于有沙波存在的槽道，參考高度為a=2/3Δ，其中Δ為沙波的高度.van Rijn[27]認(rèn)為對(duì)于沒有沙波的槽道，參考高度為a=ks，對(duì)于有沙波存在的槽道，參考高度為a=Δ，并在他的研究中a≥0.01H，H為水深.在本文中采用van Rijn提出的方法確定參考高度.在LES中，由于在近床面處的網(wǎng)格較小，因此懸沙濃度從離開床面的第n個(gè)網(wǎng)格開始計(jì)算，如圖1中加粗的網(wǎng)格所示.

圖1 計(jì)算網(wǎng)格示意圖Fig.1 The sketch map of the computational grids

3 計(jì)算算例和邊界條件

在本文中設(shè)置 2個(gè)算例，雷諾數(shù)Re均為20000，建立相應(yīng)于流向、展向和垂向的x×y×z的坐標(biāo)系.算例1為循環(huán)算例，計(jì)算區(qū)域大小為6H×3H×H(H表示水深)，x，y和z向的網(wǎng)格數(shù)分別為210，210和126.x和y向采用均勻網(wǎng)格，采用摩阻流速無量綱化后網(wǎng)格尺寸分別為Δx+=28.08和Δy+=14.04，z向采用不均勻網(wǎng)格，從逐漸擴(kuò)大到在計(jì)算中首先計(jì)算120個(gè)無量綱時(shí)間(H/U)，待計(jì)算穩(wěn)定后，統(tǒng)計(jì)了300個(gè)無量綱時(shí)間，在大型并行計(jì)算機(jī)上采用36個(gè)核計(jì)算完成.

算例1中，在流向上和展向上采用循環(huán)邊界條件.水流計(jì)算中在自由水面處引入剛蓋假定，施加滑移邊界條件.泥沙離開或進(jìn)入自由水面的凈通量為零，在自由水面處采用零通量邊界條件 (式 (6)).在床面處水流計(jì)算中施加無滑移條件，采用Werner等[28]提出的壁函數(shù)計(jì)算最靠近壁面的第一層網(wǎng)格的流速值

假定LES中近底部瞬時(shí)泥沙濃度符合與RANS相同的計(jì)算公式，同時(shí)為了使得計(jì)算不發(fā)散，引入εl有

其中sb和sn分別為為近網(wǎng)格底部和中心的泥沙濃度，sb?為近網(wǎng)格底部的挾沙力，zb和zn分別為近網(wǎng)格底部和中心的高程.εl為大尺度運(yùn)動(dòng)擴(kuò)散系數(shù)，表示大渦模擬能夠識(shí)別的包含在大尺度運(yùn)動(dòng)中的擴(kuò)散系數(shù)

其中，-u′w′為zb處的切應(yīng)力，?u/?z為相應(yīng)位置瞬時(shí)流向流速u的z向梯度.σl為大尺度運(yùn)動(dòng)Schmidt數(shù)，取為1.為了保證計(jì)算穩(wěn)定性，大尺度運(yùn)動(dòng)擴(kuò)散系數(shù)εl采用統(tǒng)計(jì)平均值

其中〈〉表示統(tǒng)計(jì)平均.

算例2為單向流動(dòng)算例，計(jì)算區(qū)域大小為48H× 3H×H，分為先導(dǎo)計(jì)算域和主計(jì)算域，其中先導(dǎo)計(jì)算域的大小為6H×3H×H，在先導(dǎo)計(jì)算域中只計(jì)算水流，不計(jì)算泥沙.在算例2中，x,y和z向的網(wǎng)格數(shù)分別為1680,210和126.網(wǎng)格尺寸同算例1.在計(jì)算中首先計(jì)算288個(gè)無量綱時(shí)間，待計(jì)算穩(wěn)定后，統(tǒng)計(jì)了1800個(gè)無量綱時(shí)間，在大型并行計(jì)算機(jī)上采用288個(gè)核計(jì)算完成.

算例2中，水流計(jì)算邊界條件的設(shè)置與算例1相同.在主計(jì)算域中計(jì)算水流和泥沙，展向和垂向邊界條件與算例1相同.流向上的邊界條件設(shè)置為：在主計(jì)算域的進(jìn)口處，設(shè)置進(jìn)口邊界條件，使得在先導(dǎo)計(jì)算域中充分發(fā)展的湍流順利進(jìn)入主計(jì)算域；在出口處，泥沙和水流都采用對(duì)流邊界條件，并且實(shí)時(shí)校正水流流量，保證流量守恒.詳見圖2.

圖2 算例2計(jì)算區(qū)域平面示意圖Fig.2 The sketch map of the computational domain in case 2

4 模型的率定與驗(yàn)證

利用算例1的水流計(jì)算結(jié)果對(duì)LES模型進(jìn)行率定和驗(yàn)證.采用算例 1的結(jié)果和直接數(shù)值模擬(direct numerical simulation,DNS)[29]計(jì)算結(jié)果進(jìn)行比對(duì)，流向流速和湍動(dòng)強(qiáng)度與DNS的比較分別見圖3(a)和圖3(b).在圖3(a)中縱坐標(biāo)為u+，即對(duì)流向的平均流速采用摩阻流速u?進(jìn)行無量綱化處理有u+=〈u〉/u?，〈u〉表示對(duì)流向流速取統(tǒng)計(jì)平均，摩阻流速u?=0.04915U，U為斷面平均流速.橫坐標(biāo)為z+，采用摩阻流速u?和υ進(jìn)行無量綱化.在計(jì)算中，基于摩阻流速u?的雷諾數(shù)為Re?=u?H/υ=983，DNS中Re?=950.從圖中可以看出，LES計(jì)算得到的結(jié)果與DNS結(jié)果符合良好，僅在自由水面處有細(xì)微的差別，LES計(jì)算得到的u+比DNS結(jié)果偏小0.8.不難發(fā)現(xiàn)，當(dāng)z+介于[30，200]時(shí)，LES計(jì)算得到的u+服從對(duì)數(shù)律分布，這與以往的研究結(jié)果是一致的.

圖3 算例1中LES和DNS流向速度和湍動(dòng)強(qiáng)度比較Fig.3 Comparisons of the stream wise velocity and turbulence intensities between LES and DNS in case 1

圖3 算例1中LES和DNS流向速度和湍動(dòng)強(qiáng)度比較(續(xù))Fig.3 Comparisons of the stream wise velocity and turbulence intensities between LES and DNS in case 1(continued)

在圖 3(b)中比較了流向、展向和垂向上的無量綱湍動(dòng)強(qiáng)度的分布.其中從圖中可以發(fā)現(xiàn)，在整體上LES計(jì)算的湍動(dòng)強(qiáng)度與DNS結(jié)果符合較好，在床面附近，LES計(jì)算得到的的最大值為2.96，比DNS偏大0.14.和w′+的最大值分別為1.33和1.04，比DNS分別偏小0.12和0.04.在z+介于[200,800]時(shí)，LES計(jì)算得到的略微偏小，而與DNS結(jié)果基本一致.由于LES中在自由水面處采用剛蓋假定，并施加滑移(對(duì)稱)邊界條件，因此水面處的而略微上抬.

圖4 無量綱切應(yīng)力計(jì)算結(jié)果對(duì)比圖Fig.4 The comparison of the dimensionless shear-stresses between LES and DNS

5 結(jié)果分析

5.1 時(shí)均的水流結(jié)果分析

由于循環(huán)槽道的流速分布等已有了詳細(xì)的討論，這里集中分析長(zhǎng)槽道算例(算例2)中的水流結(jié)果.由于在長(zhǎng)槽道的展向y上采用循環(huán)邊界條件，水流在展向上是均勻的，因此分析采用展向平均的結(jié)果.

圖5(a)中給出的是長(zhǎng)槽道中x=0,10H,20H, 30H,40H處的瞬時(shí)流向流速u/U和時(shí)均流向流速〈u〉/U沿著水深的分布情況，均做展向平均處理，其中細(xì)實(shí)線表示瞬時(shí)流向流速u/U，粗實(shí)線表示時(shí)均流向流速〈u〉/U.從圖中可以看出，不同位置處時(shí)均流向流速〈u〉/U相同，瞬時(shí)流向流速u/U沿x軸也差別不大.時(shí)均流向流速〈u〉/U沿x方向均勻分布，在水面處取得最大值1.15，采用摩阻流速u?無量綱化后為23.3，與圖3(a)的結(jié)果一致.

圖5(b)中給出的是長(zhǎng)槽道中x=0,10H,20H, 30H,40H處展向平均的湍動(dòng)強(qiáng)度分布情況，其中從圖中可以看出沿著流向各個(gè)位置的湍動(dòng)強(qiáng)度基本重合，湍動(dòng)強(qiáng)度在床面處為0，隨著水深急劇增加，在近底部取得最大值3.0，1.32和1.04，之后隨著水深的增加不斷減小.在水面處，垂向湍動(dòng)強(qiáng)度為0，相應(yīng)地，流向和展向湍動(dòng)強(qiáng)度略有抬升，在總體上，長(zhǎng)槽道中湍動(dòng)強(qiáng)度的分布與算例1中的湍動(dòng)強(qiáng)度分布基本相同(圖3(b)).在算例2中，長(zhǎng)槽道單向流為均勻流動(dòng).

圖5 算例2中不同位置上流向速度和湍動(dòng)強(qiáng)度分布Fig.5 The profile of the stream-wise velocities and turbulence intensities at dif f erent locations in case 2

圖6中給出的是長(zhǎng)槽道中x=0,10H,20H,30H, 40H處的無量綱切應(yīng)力分布，從圖中可以看出，沿程的剪切應(yīng)力基本一致.在底面和水面處均為0，離開床面處，迅速減小，-0.89為其最小值，之后沿著水深基本上呈現(xiàn)線性增加的趨勢(shì)，直到水面處歸零，與算例1中結(jié)果基本相同.

5.2 時(shí)均泥沙濃度計(jì)算結(jié)果分析

圖7(a)中給出的是循環(huán)槽道中 (算例 1)無量綱的時(shí)均泥沙濃度沿著水深的分布，其中橫坐標(biāo)為〈s〉/sb，sb為輸沙平衡時(shí)近底部的參考泥沙濃度.算例1可以給出輸沙平衡時(shí)泥沙濃度沿著水深方向的分布情況.從圖中可以看出，在輸沙平衡時(shí)LES計(jì)算得到的泥沙濃度符合Rouse公式的分布.

圖6 算例2中無量綱切應(yīng)力的分布Fig.6 The profile of the dimensionless shear stresses at dif f erent locations in case 2

圖7 時(shí)均泥沙濃度分布Fig.7 The profile of time-averaged sediment concentration

在圖7(b)中，從左到右依次為：算例2長(zhǎng)槽道中x=0,10H,20H,30H,40H處的展向平均的〈s〉/sb和算例1中輸沙平衡時(shí)統(tǒng)計(jì)得到的〈s〉/sb.從圖中可以看出在長(zhǎng)槽道中，泥沙濃度始終小于輸沙平衡時(shí)的泥沙濃度〈s〉/sb，可以判斷出長(zhǎng)槽道處于凈沖刷狀態(tài).沿著x方向，隨著沖刷的進(jìn)行，泥沙濃度由于受到水流切應(yīng)力的作用，泥沙顆粒不斷起動(dòng)，在湍動(dòng)擴(kuò)散的作用下，泥沙在不斷向前運(yùn)動(dòng)的同時(shí)，也向上擴(kuò)散，〈s〉/sb沿程逐漸增大.

圖8中給出了算例2長(zhǎng)槽道中瞬時(shí)和時(shí)均的斷面平均泥沙濃度的沿程分布，從圖中可以看出，雖然經(jīng)過斷面平均，瞬時(shí)的泥沙濃度仍然有較大的脈動(dòng).但在總體趨勢(shì)上，瞬時(shí)的斷面平均泥沙濃度呈現(xiàn)與時(shí)均的斷面平均的泥沙濃度一致的趨勢(shì)，泥沙濃度隨著x的增大而逐漸變大.從圖中可以看出，從x=0H開始，斷面平均泥沙濃度〈s〉/sb沿著x方向逐漸增大，但增大速度逐漸減慢，在x=42H處，達(dá)到最大值0.14，但由于斷面平均泥沙濃度〈s〉/sb關(guān)于x的導(dǎo)數(shù)大于零，整個(gè)長(zhǎng)槽道中，床面一直處于凈沖刷狀態(tài).

圖8 算例2中斷面平均的泥沙濃度沿程分布Fig.8 Distribution of the cross-section averaged sediment concentration in case 2

泥沙濃度湍動(dòng)強(qiáng)度沿水深的分布見圖9(a).圖中橫坐標(biāo)z采用摩阻流速u?和黏性系數(shù)υ進(jìn)行無量綱化，有z+=zu?/υ，采用平衡時(shí)近底部的泥沙濃度sb將泥沙湍動(dòng)強(qiáng)度s′進(jìn)行無量綱化.從圖中可以看出泥沙的湍動(dòng)強(qiáng)度s′/sb在近底面處達(dá)到最大值0.49，之后隨著水深的增大而急劇降低，在水面附近僅為0.04，在自由水面處幾乎為零.圖9(b)給出了采用LES得到的垂向上泥沙的湍動(dòng)通量沿水深的分布.橫坐標(biāo)為z+=zu?/υ，縱坐標(biāo)w′s′采用u?sb進(jìn)行無量綱化，在近底部w′s′/u?sb取得最大值0.115，之后隨著水深的增大不斷減小，在水面處基本為零.

圖9 算例1中泥沙濃度湍動(dòng)強(qiáng)度和湍動(dòng)通量的分布Fig.9 The turbulence intensity of sediment concentration and vertical sediment turbulence flu in case 1

5.3 湍動(dòng)黏性系數(shù)與擴(kuò)散系數(shù)分布

對(duì)于湍動(dòng)切應(yīng)力和通量有

其中vt，εt分別為湍動(dòng)黏性系數(shù)和湍動(dòng)擴(kuò)散系數(shù)，是湍動(dòng)切應(yīng)力是垂向湍動(dòng)泥沙通量.

圖10中給出了循環(huán)槽道中(算例1)和長(zhǎng)槽道中(算例2)采用LES得到的湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)沿著水深的分布.圖中的橫坐標(biāo)z采用水深H進(jìn)行無量綱化，縱坐標(biāo)湍動(dòng)黏性系數(shù)υt和湍動(dòng)擴(kuò)散系數(shù) εt，分別采用分子擴(kuò)散系數(shù) υ和 ε進(jìn)行無量綱化，υ=1/20000，取分子Schmidt數(shù)為1，因此ε=1/20000.湍動(dòng)擴(kuò)散系數(shù)εt理論值采用式(9)計(jì)算

其中κ為卡帕常數(shù)，u?為摩阻流速，z為離開床面的距離，H為水深.

從圖10(a)中可以看出，采用大渦模擬計(jì)算得到的湍動(dòng)黏性系數(shù)υt/υ關(guān)于z/H=0.5不對(duì)稱分布，在z/H=0.35處，υt/υ取得最大值60，而大渦模擬計(jì)算得到的擴(kuò)散系數(shù)εt/ε關(guān)于z/H=0.5對(duì)稱分布，在水深中間處，εt/ε取得最大值70.在z/H介于[0.15, 0.85]的區(qū)間中時(shí)，εt/ε的理論值明顯大于LES計(jì)算值.在水深中間處由理論公式得到的εt/ε的最大值是LES計(jì)算得到的εt/ε最大值的1.4倍.這是因?yàn)槔碚摴酵茖?dǎo)時(shí)假定：(1)流速服從對(duì)數(shù)分布；(2)泥沙湍動(dòng)擴(kuò)散系數(shù)等于動(dòng)量擴(kuò)散系數(shù).而實(shí)際上流速只在對(duì)數(shù)區(qū)服從對(duì)數(shù)分布，由于水流中泥沙的跟隨性不如普通的污染物，計(jì)算得到的湍動(dòng)擴(kuò)散系數(shù)與理論值有差別.

圖10 湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)分布圖Fig.10 The turbulence viscosity coefficient and the turbulence dif f usion coefficient in two cases

圖10(b)給出了長(zhǎng)槽道中(算例2)x=10H,20H, 30H,40H處展向平均的湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)的分布.從圖中可以看出，由于統(tǒng)計(jì)樣本數(shù)較少，因此動(dòng)黏性系數(shù)υt和擴(kuò)散系數(shù)εt有輕微的波動(dòng)，但總體上看，黏性系數(shù)υt和擴(kuò)散系數(shù)εt沿著x向的各個(gè)斷面基本相同，并且沿著垂向呈現(xiàn)先變大后變小的趨勢(shì).

圖11 算例1中湍動(dòng)Schmidt數(shù)分布圖Fig.11 The profil of the turbulence Schmidt number in case 1

圖11給出了湍動(dòng) Schmidt數(shù) σt沿著水深的分布.從圖中可以看出σt沿著水深不斷變化，湍動(dòng)Schmidt數(shù)σt在接近水面和底部時(shí)較大，在水深中間處較小，當(dāng)z/H介于[0.25,0.85]中時(shí)，σt小于1.在z/H=0.55時(shí)，σt=0.83為最小值.但在雷諾時(shí)均模擬中σt常常假定為常數(shù)[30-31]，當(dāng)水流結(jié)果計(jì)算準(zhǔn)確時(shí)，σt存在偏差，從而在泥沙濃度計(jì)算中引入誤差，尤其對(duì)近底部和水面附近的泥沙濃度有較大的影響.

6 結(jié)論

在本文中采用基于并行計(jì)算的大渦模擬程序，借鑒了雷諾時(shí)均模擬中懸沙計(jì)算中的參考高度概念，建立了LES中的細(xì)顆粒泥沙凈沖刷條件下的近底邊界條件模擬了循環(huán)槽道和長(zhǎng)槽道中水流運(yùn)動(dòng)和泥沙輸移的過程.在輸沙平衡條件下，泥沙的濃度分布服從Rouse公式的分布.長(zhǎng)槽道流動(dòng)中，流動(dòng)沿著x向均勻分布，整個(gè)槽道處于凈沖刷狀態(tài)，在水流切應(yīng)力和湍動(dòng)作用下，泥沙顆粒不斷起動(dòng)并向上擴(kuò)散，泥沙濃度沿程不斷增大.泥沙濃度的湍動(dòng)強(qiáng)度和垂向湍動(dòng)通量在近底處達(dá)到最大，沿著垂向迅速衰減.通過LES統(tǒng)計(jì)得到的湍動(dòng)黏性系數(shù)和擴(kuò)散系數(shù)沿著水深方向呈現(xiàn)先增加后減小的趨勢(shì)，整體上湍動(dòng)擴(kuò)散系數(shù)的計(jì)算值小于實(shí)際值.湍動(dòng)Schmidt數(shù)不是常數(shù)，當(dāng)z/H介于[0.25，0.85]中時(shí)，σt小于1.在z/H=0.55時(shí)，σt=0.83為最小值.

1 Walling DE,Owens PN,Carton J,et al.Storage of sedimentassociated nutrients and contaminants in river channel and flood plain systems.Applied Geochemistry,2003,18(1)：195-220

2 Collins AL,Walling DE,Leeks GJ.Storage of fine-grainesediment and associated contaminants within the channels of lowland permeable catchments in the UK//Sediment Budgets 1,Walling D E,Horowitz A J,eds.IAHS Publication No.291,IAHS Press：Wallingford,UK,2005：259-268

3 Fang HW,Zhao HM,Shang QQ,et al.Ef f ect of biofil on the rheological properties of cohesive sediment.Hydrobiologia,2012, 694(1)：171-181

4 Huang L,Fang H,Reible D.Mathematical model for interactions and transport of phosphorus and sediment in the Three Gorges Reservoir.Water Research,2015,85：393-403

5 Fang HW,Huang L,Wang JY,et al.Environmental assessment of heavy metal transport and transformation in the Hangzhou Bay, China.Journal of Hazardous Materials,2015,302(17)：3836-3849

6 Marchioli C,Soldati A.Mechanisms for particle transfer and segregation in a turbulent boundary layer.Journal of Fluid Mechanics, 2002,468：283-315

7 Chang YS,Scotti A.Entrainment and suspension of sediments into a turbulent fl w over ripples.Journal of turbulence,2003,4(19)：1-22

8 van Rijn LC.2007.Unifieview of sediment transport by currents and waves II：suspended transport.Journal of Hydraulic Engineering-ASCE,2007,133(6)：668-689

9 Zeng J,Constantinescu G,Weber L.A 3D non-hydrostatic model to predict fl w and sediment transport in loose-bed channel bends.Journal of Hydraulic Research,2008：46(3)：356-372

10 Fang H,He G,Liu J,et al.3D numerical investigation of distorted scale in hydraulic physical model experiments.Journal of Coastal Research,2015,10052(1)：41-54

11 Keylock CJ,Hardy RJ,Parsons DR,et al.The theoretical foundations and potential for large-eddy simulation(LES)in fluvia geomorphic and sedimentological research.Earth-Science Reviews, 2005,71(3-4)：271-304

12 Rodi W.Comparison of LES and RANS calculations of the fl w around bluf fbodies.Journal of Wind Engineering and Industrial Aerodynamics,1997,69-71：55-75

13 Zedler EA,Street RL.Large-eddy simulation of sediment transport：currents over ripples.Journal of Hydraulic Engineering-ASCE, 2001,127(6)：444-452

14 Zedler EA,Street RL.Sediment transport over ripples in oscillatory fl w.Journal of Hydraulic Engineering-ASCE,2006,132(2)：180-193

15 Chou YJ,Fringer OB.Modeling dilute sediment suspension using large-eddy simulation with a dynamic mixed model.Physics of Fluids,2008,20(11)：115103

16 Widera P,Ghorbaniasl G,Lacor C.Study of the sediment transport over fla and wavy bottom using large-eddy simulation.Journal of Turbulence,2009,10(33)：1-20

17 Widera P,Toorman E,Lacor C.Large eddy simulation of sediment transport in open-channel fl w.Journal of Hydraulic Research, 2009,47(3)：291-298

18 Bai J,Fang HW,Stoesser T.Transport and deposition of fin sediment in open channels with dif f erent aspect ratios.Earth Surface Processes and Landforms,2013,38(6)：591-600

19 Germano M,Piomelli U,Moin P,et al.A dynamic subgrid-scale eddy viscosity model.Physics of Fluids,1991,3(7)：1760-1765

20 Breuer M,Rodi W.Large-eddy simulation of turbulent fl w through a straight square duct and a 180°bend//Voke PR,Kleiser L,Chollet JP,eds.Direct and Large-Eddy Simulation I.Netherlands：Kluwer, 1994.273-285

21 Hinterberger C,Frohlich J,Rodi W.Three-dimensional and depthaveraged large-eddy simulations of some shallow water fl ws.Journal of Hydraulic Engineering-ASCE,2007,133(8)：857-872

22 Stoesser T,Nikora V.Flow structure over square bars at intermediate submergence：large eddy simulation(LES)study of bar spacing ef f ect.Acta Geophysica,2008,56(3)：876-893

23 Stoesser T,Ruether N,Olsen NRB.Calculation of primary and secondary fl w and boundary shear stresses in a meandering channel.Advances in Water Resources,2010,33(2)：158-170

24 Zhu J.A low-dif f usive and oscillation-free convection scheme.Communications in Applied Numerical Methods,1991,7(3)：225-232

25 Stone HL.Iterative solution of implicit approximations of multidimensional partial dif f erential equations.SIAM Journal on Numerical Analysis,1968,5(3)：530-558

26 CelikI,RodiW.Modelingsuspendedsedimenttransportinnonequilibrium situations.Journal of Hydraulic Engineering-ASCE,1988, 114(10)：1157-1191

27 van Rijn LC.Sediment transport,part II：suspended load transport.Journal of Hydraulic Engineering-ASCE,1984,110(11)：1613-1641

28 Werner H,Wengle H.Large-eddy simulation of turbulent fl w over and around a cube in a plate channel//Proceeding of 8th Symposium on Turbulent Shear Flows,1991.155-168

29 Mito Y,Hanratty TJ,Zandonade P,et al.Flow visualization of super bursts and of the Log-Layer in a DNS at Reτ=950.Flow Turbulence and Combustion,2007,79(2)：175-189

30 Wu W,Rodi W,Wenka T.3D numerical modeling of fl w and sediment transport in open channels.Journal of Hydraulic Engineering-ASCE,2000,126(1)：4-15

31 Fang HW,Wang GQ.Three-dimensional mathematical model of suspended-sediment transport.Journal of Hydraulic Engineering-ASCE,2000,126(8)：578-592

NUMERICAL SIMULATION OF EROSION AND TRANSPORT OF FINE SEDIMENTS BY LARGE EDDY SIMULATION

Bai Jing?Fang Hongwei?,1)He Guojian?Xie Chongbao?Gao Hong??

?(China Irrigation and Drainage Development Center，Beijing100054，China)

?(State Key Laboratory of Hydro Science and Engineering,Department of Hydraulic Engineering，Tsinghua University，Beijing100084，China)

??(China Green Water International Consulting Co.，Ltd，Beijing 100054，China)

In general Reynolds-averaged simulation(RANS)is used in the traditional numerical simulation of water fl w and sediment transport.Large eddy simulation(LES)can reflec fl w structures more accurately and give more details of water fl w compared with RANS.The development of computing technology makes it possible to study the rules of water fl w and sediment transport by an LES model.In this paper,we tried to introduce boundary conditions for suspended sediment transport for the LES model under the net-erosion condition.Water fl w and sediment transport in a cyclic case and a one-way fl w case were calculated via the LES model with a dynamic sub-grid stresses module and a suspended sediment calculation module in the paper.Direct numerical simulation(DNS)results were used to calibrate the LES model and the results from LES showed good agreements with the DNS results.The distribution characteristics of sediment concentration,turbulence intensity and turbulent flu es of sediment were explored in the paper.Under the net-erosion condition,the equilibrium sediment concentration profil was coincident with the line of the Rouse equation.It showed that the turbulence intensity and turbulent flu es of sediment had peak values near the bottom and then decreased rapidly along the vertical direction.The turbulent viscosity and dif f usion coefficients were calculated and their peak values were or near the mid-depth of water.The turbulent Schmidt number was not constant along the vertical direction,and it was larger near the free surface and the bottom while it was smaller near the mid-depth of water fl w.

fin sediment,net erosion,LES,sediment transport

TV142

A doi：10.6052/0459-1879-16-235

2016-08-25收稿，2016-11-15錄用,2016-11-18網(wǎng)絡(luò)版發(fā)表.

1)方紅衛(wèi)，教授，主要研究方向：環(huán)境泥沙及河流動(dòng)力學(xué).E-mail：fanghw@mail.tsinghua.com

白靜，方紅衛(wèi)，何國(guó)建，謝崇寶，高虹.細(xì)顆粒泥沙凈沖刷和輸移的大渦模擬研究.力學(xué)學(xué)報(bào)，2017,49(1)：65-74

Bai Jing,Fang Hongwei,He Guojian,Xie Chongbao,Gao Hong.Numerical simulation of erosion and transport of fin sediments by large eddy simulation.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：65-74