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

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

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

?(清華大學水利水電工程系，水沙科學與水利水電工程國家重點實驗室，北京100084)

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

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

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

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

?(清華大學水利水電工程系，水沙科學與水利水電工程國家重點實驗室，北京100084)

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

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

細顆粒泥沙,凈沖刷,大渦模擬,泥沙輸移

引言

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

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

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

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

1 控制方程

水流計算的控制方程為

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

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

懸沙運動的控制方程為

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

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

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

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

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

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

3 計算算例和邊界條件

在本文中設(shè)置 2個算例，雷諾數(shù)Re均為20000，建立相應(yīng)于流向、展向和垂向的x×y×z的坐標系.算例1為循環(huán)算例，計算區(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)格，從逐漸擴大到在計算中首先計算120個無量綱時間(H/U)，待計算穩(wěn)定后，統(tǒng)計了300個無量綱時間，在大型并行計算機上采用36個核計算完成.

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

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

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

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

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

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

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

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

4 模型的率定與驗證

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

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

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

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

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

5 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

白靜，方紅衛(wèi)，何國建，謝崇寶，高虹.細顆粒泥沙凈沖刷和輸移的大渦模擬研究.力學學報，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