后臺階下游水流泥沙運動模擬研究

孫建軍，張慶河

（天津大學(xué)水利工程仿真與安全國家重點實驗室，天津300072）

海岸河口及港口工程

后臺階下游水流泥沙運動模擬研究

孫建軍，張慶河

（天津大學(xué)水利工程仿真與安全國家重點實驗室，天津300072）

利用計算流體力學(xué)軟件OpenFOAM和離散元軟件LIGGGHTS建立的水流-泥沙運動耦合模型對后臺階下游三維水流流動及泥沙的起動和輸運進行了數(shù)值模擬。結(jié)果表明，動態(tài)一方程紊流模型能較好模擬臺階下游時均流場、紊動強度以及雷諾應(yīng)力的分布，通過耦合模型計算得到的臺階下游泥沙起動概率以及床面輸沙率分布與實驗結(jié)果符合較好。計算結(jié)果揭示在臺階下游復(fù)雜水流條件下，紊流結(jié)構(gòu)及渦旋運動對泥沙起動有重要作用，而時均流動對床面輸沙起關(guān)鍵作用。

LES?DEM耦合模型；后向臺階流；動態(tài)一方程紊流模型；泥沙運動

復(fù)雜流動下的泥沙起動和輸運是河流、海岸工程研究的重點課題，一直受到人們的關(guān)注。目前，復(fù)雜流動下的泥沙輸運問題通常依賴于經(jīng)驗公式，往往通過床面切應(yīng)力公式或挾沙力公式得到輸沙率［1-3］來進行計算。近年來，隨著離散元法（DEM）的發(fā)展，運用計算流體力學(xué)（CFD）與離散元法（DEM）耦合數(shù)值模型，從細觀尺度模擬水流泥沙運動得到了快速發(fā)展［4-6］，能夠刻畫大量泥沙顆粒在近床面運動的細節(jié)，它有可能揭示復(fù)雜流動下泥沙顆粒的運動規(guī)律。蘇東升等［7］基于計算流體力學(xué)開源軟件OpenFOAM、離散顆粒開源軟件LIGGGHTS以及CFDEM耦合庫［8-9］，建立了描述水流泥沙運動的離散顆粒模型，從細觀尺度研究了明渠流中近床面泥沙的一般運動規(guī)律，但將其應(yīng)用于復(fù)雜流動下泥沙運動規(guī)律研究的工作幾乎沒有，本文將在上述工作的基礎(chǔ)上研究后臺階下游的三維水流與泥沙運動規(guī)律。

1 LES?DEM耦合模型

1.1流體控制方程及紊流模型

1.1.1 控制方程

OpenFOAM中，考慮顆粒組分影響的三維不可壓流體運動的控制方程由以下連續(xù)性方程和動量方程組成

式中：αf為流體所占據(jù)的體積分數(shù)；f為流體的速度向量；p為流體壓力；為重力加速度向量；為流體有效應(yīng)力的張量；pf為顆粒與流體之間交換的動量項。

1.1.2 紊流模型

為了使整個控制方程封閉，運用大渦模擬（LES）方法，LES方程向N?S動量方程中引入了額外的應(yīng)力，稱為亞格子尺度應(yīng)力（簡稱SGS應(yīng)力），為求解SGS應(yīng)力，其中以渦粘假設(shè)為基礎(chǔ)的理論公式如下

Samgorinsky模型是目前廣泛應(yīng)用的亞格子模型，但由于其采用單一常數(shù)Cs計算νsgs，會存在耗散過大問題［10］。因此本文采用動態(tài)一方程模型，即計算系數(shù)隨空間和時間調(diào)整的動態(tài)模型

式中：ksgs代表亞格子紊動動能，采用以下控制方程求解

1.2 顆粒運動控制方程

LIGGGHTS顆粒模型中，離散顆粒的運動遵循牛頓第二定律，其控制方程包含平動和轉(zhuǎn)動兩部分［11］

為準確合理地描述顆粒運動，確定顆粒之間和顆粒與固壁間受力以及顆粒所受的流體作用力是模型的關(guān)鍵。蘇東升［7］的研究結(jié)果表明，顆粒之間和顆粒與固壁間的接觸碰撞采用Hertz?Mindlin無滑移軟球模型[12]，流體作用于顆粒上的拖曳力采用Benyahia拖曳力公式[13]可以較好描述上述相互作用力。

圖1 計算域（單位：m）Fig.1 Computational domain

1.3顆粒與流體運動的耦合

LES?DEM耦合模型，流體運動采用OpenFOAM中基于PISO算法的瞬態(tài)不可壓縮流體求解器求解，顆粒運動采用LIGGGHTS計算，利用基于Open?FOAM程序框架的CFDEM耦合庫完成流體和顆粒之間的相互作用力計算及信息傳遞。其計算過程分以下步驟［14］

（1）LIGGGHTS程序計算顆粒運動；

（2）CFDEM耦合庫計算顆粒所在網(wǎng)格單元的體積分數(shù)；

（3）CFDEM耦合庫計算流體與顆粒之間相互作用力，從而得到動量交換項；

（4）耦合庫計算得到的相互作用力傳遞至LIGGGHTS程序；

圖2 臺階下游平均流速剖面（Z=2.5D處）Fig.2 Mean streamwise velocity profiles at the downstream of the stepZ=2.5D

（5）基于動量交換項，OpenFOAM程序求解流體運動，并傳遞至CFDEM耦合庫；

（6）重復(fù)步驟1～5。

圖3 臺階下游時均流速矢量場圖（Z=2.5D處）Fig.3 Mean streamwise velocity vectors at the downstream of the stepZ=2.5D

2 后臺階流場模擬

為保證模型流場計算結(jié)果的合理性，首先對臺階下游床面無泥沙顆粒的情況進行數(shù)值模擬并和實驗結(jié)果進行比較以驗證模型的合理性。計算域范圍為-5D≤X≤20D，0≤Y≤5D 0≤Z≤5D，D為臺階高度，計算域如圖1所示。

入流邊界水流流速根據(jù)對數(shù)流速剖面原理給定［15］

式中：u為不同水深處的流速；u*為摩阻流速；ν表示運動粘滯系數(shù)；κ表示卡門常數(shù)，取0.41，保證模擬的入流平均流速U0=0.2 m/s，流動雷諾數(shù)ReD=U0D/ν=5 000，上邊界采用滑移邊界條件，底邊界采用無滑移邊界條件，Z軸方向設(shè)為循環(huán)邊界。上述計算條件與Jovic和Driver實驗條件一致[16]。

根據(jù)Le等[17]模擬研究，本文模擬時間取ttotal=400D/U0=50 s，流場時均統(tǒng)計時間取ΔTave=109D/U0≈14 s，可以獲得合理時均流場與紊動特征。

2.1 時均流速剖面驗證

根據(jù)大渦模擬的瞬時流場進行時均統(tǒng)計，得到臺階下游平均流速剖面。圖2給出了本文計算得到的平均流速剖面與Jovic和Driver實驗[16]的對比結(jié)果，可以看出，數(shù)值模擬結(jié)果與實驗結(jié)果一致。圖3為臺階下游時均流速矢量場圖，由圖可知，在臺階后流動會出現(xiàn)分離，形成回流區(qū)，超過這一位置，分離流重新附著床面（再附），逐漸恢復(fù)為對數(shù)剖面，可以看出X=4D位置處在回流區(qū)內(nèi)，X=6D位置在再附點附近。

2.2 紊動強度與雷諾應(yīng)力驗證

圖4～圖6顯示了后臺階下游紊動強度與雷諾應(yīng)力的計算結(jié)果與實驗結(jié)果，圖7為臺階下游紊動強度分量與雷諾應(yīng)力的三維剖面圖。由圖可知，在0.5D高度臺階位置，流向和垂向紊動強度都較強，在床面附近，流向紊動強度較大，而垂向紊動強度則較小。

圖4 臺階下游流向紊動強度剖面Fig.4 Relative streamwise turbulence intensity profiles at the downstream of the step

3 臺階下游泥沙輸運模擬

3.1 模擬參數(shù)設(shè)定

劉春嶸［18］采用圖像技術(shù)研究了后臺階下游床面泥沙的起動概率，并建立數(shù)學(xué)模型進行了數(shù)值模擬，為方便對比驗證模型，本文采用與劉春嶸實驗相同的模型布置，在臺階下游20D長度范圍內(nèi)隨機鋪設(shè)100 000個泥沙顆粒，顆粒直徑d=0.425 mm，密度 ρs=1 400 kg/m3，計算流體密度ρf=1 000 kg/m3，動 力 粘 滯 系 數(shù)μ=10-6kg/m·s，顆粒計算的Hertz?Mindlin無滑移“軟球”模型參數(shù)如表1所示，流體與顆粒計算的時間間隔分別為Δtf=10-3s和Δtp=10-5s，顆粒每計算100步與流體進行一次耦合。

本文根據(jù)泥沙運動速度的大小來判斷泥沙是否起動，在某時刻t，沿臺階下游以0.4D寬為間隔分區(qū)域統(tǒng)計不同位置處的泥沙起動概率Pit，起動的顆粒數(shù)目與該區(qū)域顆?？倲?shù)之比即為該區(qū)域t時刻的泥沙起動概率Pit，統(tǒng)計時間間隔為0.05 s，共取200個時刻值進行平均，從而得到該位置處的泥沙起動概率Pi。

圖5 臺階下游垂向紊動強度剖面Fig.5 Relative vertical turbulence intensity profiles at the downstream of the step

圖6 臺階下游雷諾應(yīng)力剖面Fig.6 Relative Reynolds stress profiles at the downstream of the step

圖8 給出了本文計算得到的臺階下游不同位置處泥沙的起動概率Pi，并與劉春嶸的實驗結(jié)果及數(shù)值模擬結(jié)果進行了比較，可以看出，本文計算結(jié)果與實驗值基本一致，在靠近臺階區(qū)域泥沙起動概率較低，隨著距臺階距離的增大，泥沙起動概率逐漸增大，在再附點振蕩區(qū)域內(nèi)達到最大值。

3.2 泥沙輸運

目前關(guān)于床面泥沙輸運問題的研究，大部分是基于平直槽道流動下的結(jié)果，一般都以床面切應(yīng)力作為泥沙起動輸運的水動力因素［19］，而在復(fù)雜流動條件下，由于渦旋以及強紊動作用的存在，只考慮平均床面切應(yīng)力因素并不合適，Nelson等[20]利用LDV測速技術(shù)對后臺階下游泥沙輸運與紊動相互作用問題進行了研究，Dou等［21］將渦量強度、紊動強度等因素等權(quán)重線性疊加，建立了輸沙率模型，這些工作都是對復(fù)雜流動下泥沙運動水動力機制的初步探索，這里采用大渦模擬準確模擬流場，并與離散顆粒模型耦合，從細觀尺度上研究復(fù)雜流動條件下泥沙的運動規(guī)律。

圖7 臺階下游紊動強度及雷諾應(yīng)力分布圖Fig.7 Turbulence intensity and Reynolds stress profiles at the downstream of the step

本文根據(jù)Nelson等人[20]的實驗布置進行算例設(shè)計，臺階高度D=0.04 m，泥沙顆粒粒徑服從高斯分布，中值粒徑d50=0.9 mm，標準差為σ=0.1 mm，密度ρs=2 650 kg/m3，泥沙顆粒隨機均勻平鋪在臺階下游床面，共計415 000個顆粒。

初始計算時，固定泥沙顆粒直至流場紊流充分發(fā)展穩(wěn)定后，釋放床面泥沙顆粒開始運動，并持續(xù)模擬30 s。為保證結(jié)果的準確性，結(jié)果統(tǒng)計分析取后20 s時間內(nèi)的數(shù)據(jù)。沿著臺階下游順流方向，以0.01 m寬為間隔分區(qū)域計算不同位置處的瞬時輸沙率qsx，再進行長時間的統(tǒng)計平均得到不同位置處的平均輸沙率qˉsx。根據(jù)泥沙顆粒的瞬時速度可以計算得到泥沙的瞬時輸沙率，公式如下

表1 模型中顆粒參數(shù)Tab.1 Particle parameters used for the present simulation

無量綱輸沙率公式如下

圖8 臺階下游不同位置處的泥沙起動概率比較Fig.8 Comparison of the sediment incipient probability at the downstream of the step

式中:up,i為第i個泥沙顆粒的速度；Vp,i為第i個泥沙顆粒的體積；A為模擬區(qū)域平面的面積；ρs為顆粒密度。

圖9～圖10給出了臺階下游床面附近（Y=0.005 m）不同位置處的平均流速及流向紊動強度，圖11為臺階下游縱向時均流速矢量場圖，圖12給出了后臺階下游的平均無量綱輸沙率qˉ*，并與Nelson的實驗結(jié)果以及Schmeeckle模擬結(jié)果［22］進行了比較，可以看出，在回流區(qū)內(nèi)，平均輸沙率為負，床面泥沙向上游輸運，在平均再附點（X/D≈6）附近，雖然由于大尺度渦的存在導(dǎo)致紊動強度較大，泥沙起動概率較大，但附近平均流速接近為零，故造成床面平均輸沙率較低。

圖9 臺階下游不同位置處的平均流速（Y=0.005 m）Fig.9 Mean streamwise velocity at the downstream of the stepY=0.005 m

圖10 臺階下游不同位置處的流向紊動強度（Y=0.005 m）Fig.10 Streamwise turbulence intensity at the downstream of the stepY=0.005 m

圖11 臺階下游時均流速矢量場圖（Z=0.05 m處）Fig.11 Mean streamwise velocity vectors at the downstream of the stepZ=0.05 m

圖12 臺階下游不同位置處的平均無量綱輸沙率比較Fig.12 Comparisons of mean dimensionless sediment transport rates at the downstream of the step

4 結(jié)論

本文利用計算流體力學(xué)軟件OpenFOAM和離散元軟件LIGGGHTS建立的水流-泥沙運動耦合模型，從細觀尺度研究了后臺階下游三維紊流流動以及床面泥沙的起動和輸運內(nèi)部機制，主要得到以下結(jié)論：

（1）動態(tài)一方程紊流模型能較為準確地模擬后臺階下游復(fù)雜的三維紊流流動，計算得到的時均流速、紊動強度及雷諾應(yīng)力剖面與實驗結(jié)果符合較好，保證了模型水動力條件的合理性。

（2）分析得到0.5D高度臺階處，流向和垂向紊動強度都較強，在床面附近，流向紊動強度較大，而垂向紊動強度較小。

（3）模擬得到的臺階下游平均輸沙率與實驗結(jié)果較為一致，回流區(qū)內(nèi)平均輸沙率為負，床面泥沙向上游輸運，平均再附點附近，雖然大尺度渦的存在導(dǎo)致紊動強度較大，泥沙起動概率較大，但由于附近平均流速接近為零，造成平均輸沙率較低。

（4）在后臺階下游復(fù)雜流動條件下，渦旋及紊動作用對泥沙起動起主要作用，而對泥沙輸運起重要作用的是時均流動。

［1］Niemann S L，F(xiàn)reds?e J，Jacobsen N G.Sand dunes in steady flow at low Froude numbers：Dune height evolution and flow resis?tance［J］.Journal of Hydraulic Engineering，2010，137（1）：5-14.

［2］Chou Y J，F(xiàn)ringer O B.A model for the simulation of coupled flow?bed form evolution in turbulent flows［J］.Journal of Geophysical Research：Oceans，2010，115（C10）.

［3］Nguyen Q，Wells J C.A numerical model to study bedform development in hydraulically smooth turbulent flows［J］.Journal of Hy?draulic Engineering，2009，53（1）：157-162.

［4］王光謙，孫其誠.顆粒物質(zhì)及其多尺度結(jié)構(gòu)統(tǒng)計規(guī)律［J］.工程力學(xué)，2009，26（S2）：1-7. WANG G Q，SUN Q C.Granular matter and the scaling laws［J］.Engineering Mechanics，2009，26（S2）：1-7.

［5］Zhu H P，Zhou Z Y，Yang R Y，et al.Discrete particle simulation of particulate systems：a review of major applications and findings［J］.Chemical Engineering Science，2008，63（23）：5 728-5 770.

［6］肖鋒軍，郭烈錦，王躍社，等.三維混合沙輸運數(shù)值模擬［J］.工程熱物理學(xué)報，2012，33（2）：259-262. XIAO F J，GUO L J，WANG Y S，et al.A 3?dimensional simulation of mixed sand transport［J］.Journal of Engineering Thermophys?ics，2012，33（2）:259-262.

［7］蘇東升，張慶河，孫建軍，等.基于CFD-DEM耦合方法的近床面水流泥沙運動模擬研究［J］.水道港口，2016，37（3）：224-230. SU D S，ZHANG Q H，SUN J J，et al.Simulation of fluid?sediment particle motion near bed based on CFD?DEM coupling method［J］.Journal of Waterway and Harbor，2016，37（3）：224-230.

［8］Sadaghiani M R S，Jentsch H，F(xiàn)aulstich K，et al.DEM modeling and identification of representative element volume of soil skeleton［J］.Numerical Methods in Geotechnical Engineering，2014：403.

［9］Zhao J，Shan T.Coupled CFD?DEM simulation of fluid?particle interaction in geomechanics［J］.Powder technology，2013，239：248-258.

［10］張海濤，曹曙陽.基于動態(tài)亞格子模型的方柱繞流大渦模擬［J］.沈陽建筑大學(xué)學(xué)報:自然科學(xué)版，2013，29（3）:434-439. ZHANG H，CAO S Y.Large Eddy Simulation of Flow Past a Square Cylinder Based on Dynamic Sub?Grid Scale Models［J］.Jour?nal of Shenyang Jianzhu University，2013，29（3）:434-439.

［11］胡國明.顆粒系統(tǒng)的離散元素法分析仿真［M］.武漢：武漢理工大學(xué)出版社，2010.

［12］孫其誠，光謙.顆粒物質(zhì)力學(xué)導(dǎo)論［M］.北京：科學(xué)出版社，2009.

［13］Benyahia S，Syamlal M，O'Brien T J.Extension of Hill?Koch?Ladd drag correlation over all ranges of Reynolds number and solids volume fraction［J］.Powder Technology，2006，162（2）：166-174.

［14］Goniva C，Kloss C，Hager A，et al.An open source CFD?DEM perspective［J］.Proc Openfoam Workshop，2010.

［15］劉春晶，李丹勛，王興奎.明渠均勻流的摩阻流速及流速分布［J］.水利學(xué)報，2005，36（8）：0950-0955. LIU C J，LI D X，WANG X K.Experimental study on friction velocity and velocity profile of open channel flow［J］.Journal of Hy?draulic Engineering，2005，36（8）:0950-0955.

［16］Jovic S，Driver D M.Backward?facing step measurements at low Reynolds number，Reh=5000［J］.NASA Tech Mem，1994:94.

［17］Le H，Moin P，Kim J.Direct numerical simulation of turbulent flow over a backward?facing step［J］.Journal of Fluid Mechanics，1997，330（1）：349-374.

［18］劉春嶸.復(fù)雜流動下底床局部沖刷實驗和數(shù)值模擬研究［D］.北京:中國科學(xué)院力學(xué)研究所，2003.

［19］Rijn L C V.Sediment Pick?Up Functions［J］.Journal of Hydraulic Engineering，1984，110（10）:1 494-1 502.

［20］Nelson J M，Shreve R L，McLean S R，et al.Role of near?bed turbulence structure in bed load transport and bed form mechanics［J］.Water Resources Research，1995，31（8）：2 071-2 086.

［21］Dou X B.Using A 3?D Model to predict local scour［C］//Steven R.Water Resource Engineering 98.Memphis，Tennessee，ASCE，1998:198-203.

［22］Schmeeckle M W.The role of velocity，pressure，and bed stress fluctuations in bed load transport over bed forms：numerical simu?lation downstream of a backward?facing step［J］.Earth Surface Dynamics，2015，3（1）：105.

Simulation of fluid?sediment particle motion at the downstream of backward?facing step flow

SUN Jian?jun，ZHANG Qing?he
（State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China）

Using computational fluid dynamics software OpenFOAM and particle motion simulation software LIGGGHTS，a coupled fluid?particle model was presented and applied in the investigation of three?dimensional tur?bulence，sediment incipience and sediment transport at the downstream of the backward?facing step flow.The simu?lation results show that the dynOneEqEddy LES model can describe the mean velocity profile，the turbulence inten?sity and the distribution of Reynold stress well.The probability of sediment incipience and the sediment flux ob?tained by the coupled model at the downstream of the backward?facing step agree well with the experiment.The re?sults indicate that the turbulence and large?scale vortices induce the sediment incipience and the mean flow contrib?utes to its transport in the case of complex flows at the downstream of the backward?facing step.

LES?DEM coupled model;backward?facing step flow;dynOneEqEddy LES model;sediment mo?tion

TV 142；O 242.1

A

1005-8443（2017）02-0109-06

2016-11-03；

2017-01-13

國家自然科學(xué)基金資助（51179122）

孫建軍（1990-），男，山東省濰坊人，碩士研究生，主要從事港口海岸及近海工程研究。

Biography：SUN Jian?jun（1990-），male，master student.