基于CFD-DEM耦合方法的近床面水流泥沙運動模擬研究

蘇東升，張慶河，孫建軍，李明星

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

基于CFD-DEM耦合方法的近床面水流泥沙運動模擬研究

蘇東升，張慶河，孫建軍，李明星

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

基于CFD-DEM方法，采用流體計算軟件OpenFOAM、顆粒運動模擬軟件LIGGGHTS及CFDEM耦合庫，建立流體—顆粒運動耦合模型，并利用模型研究不同水動力條件下近床面流體與顆粒的運動規(guī)律。模擬結(jié)果表明，模型能較好模擬近床面水流紊動特性及雷諾應力分布，也能較準確刻畫出顆粒未起動、推移質(zhì)輸沙及推移質(zhì)和懸移質(zhì)共同輸沙的三種狀態(tài)，計算得到的推移質(zhì)單寬輸沙率與實驗結(jié)果及經(jīng)驗公式吻合較好。

CFD-DEM耦合模型；近床面水流；紊流強度；推移質(zhì)輸沙

水流作用下的近床面泥沙運動是河流動力學研究的重要課題，一直受到國內(nèi)外學者的重視，目前已有大量理論、實驗及數(shù)值模擬的文獻發(fā)表［1-2］。近年來，計算流體力學（CFD）的廣泛應用及離散單元法（DEM）的日臻完善，使得水流泥沙運動研究從宏觀研究進入細觀尺度的數(shù)值模擬研究［3-5］。CFD與DEM的耦合數(shù)值模型可以描述泥沙顆粒在水體中的運動過程，從細觀尺度刻畫大量泥沙顆粒在近床面運動的細節(jié)，并可以在此基礎上進一步總結(jié)近床面水流泥沙宏觀運動規(guī)律［6-7］。與直接數(shù)值模擬（DNS）［8］、浸潤邊界法（IBM）［9］及格子玻爾茲曼（LBM）［10］等可以描述顆粒運動受力的方法相比，CFD-DEM耦合模型對于顆粒本身并不作解析，減少了對計算資源的需求，因而可以應用于描述局部沖刷等顆粒數(shù)目上百萬的實際問題，具有較為廣泛的應用前景。為此，本文將基于計算流體力學開源軟件OpenFOAM?和離散顆粒開源軟件LIGGGHTS?以及二者耦合的開源軟件CFDEM?，建立描述水流泥沙運動的細觀模型，研究不同水流動力條件下近床面水體的紊動特性與泥沙顆粒的運動規(guī)律。

1 CFD-DEM耦合模型

1.1 模型簡介

CFDEM是基于Linux平臺，采用C++語言編寫的一款面向?qū)ο蟮腃FD-DEM開源耦合庫。CFDEM耦合庫將OpenFOAM模擬的流體運動與LIGGGHTS計算的顆粒運動之間以傳遞動量項的方式進行耦合，建立了流體與顆粒的模擬框架［11］。

1.2 流體運動控制方程

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

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

1.3 顆粒運動控制方程

在LIGGGHTS顆粒模型中，顆粒運動遵循牛頓第二定律，其控制方程主要包括平動與轉(zhuǎn)動兩方面，具體形式為

1.4 耦合過程

在耦合過程中，OpenFOAM采用大渦模擬（LES）來描述流場中的紊動特性，亞格子應力采用Smagorinsky模型［12］；LIGGGHTS中采用Hertz?Mindlin無滑移“軟球”模型［13］計算接觸顆粒之間的作用力。具體的耦合過程如圖1所示。

圖1 顆粒與流體運動的耦合過程Fig.1 Procedure for coupling of particles and fluid

2 近床面水流泥沙運動的數(shù)值模擬

2.1計算區(qū)域及模擬參數(shù)設定

流體的計算域為長方體，大小為x×y×z=0.12 m×0.06 m×0.04 m。流體的流動方向為x軸正方向，y軸為展向，z軸為垂向。計算網(wǎng)格為結(jié)構(gòu)化網(wǎng)格，水平和垂向網(wǎng)格步長為1.0 mm。顆粒運動的計算域與流體重合，顆粒簡化為直徑為0.5 mm的球體。流體與顆粒計算的時間間隔分別為Δtf=10-3s和Δts=10-5s，顆粒計算100步后與流體進行一次耦合。計算流體為清水，其密度ρf=1 000 kg/m3，動力粘滯系數(shù)μ=10-6kg/m·s。用于顆粒計算的Hertz?Mindlin無滑移“軟球”模型參數(shù)如表1所示。

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

2.2算例設計及模型初邊值

為了研究不同流速下流體與泥沙運動的規(guī)律，本文共設計九組算例。算例1～9中，流體沿x方向的目標平均流速＜＜u→f＞＞從0.2 m/s增加至1.0 m/s，間隔為0.1 m/s。為了保證計算域底邊界處泥沙顆粒不被全部沖起，對于目標平均流速＜＜u→f＞＞=0.2～0.6 m/s的算例，底床泥沙顆粒的厚度取為1.7 mm，共計115731個顆粒；對于目標平均流速＜＜u→f＞＞=0.7～1.0 m/s的算例，底床泥沙顆粒的厚度取為5.0 mm左右，共計330 023個顆粒。

泥沙顆粒隨機均勻平鋪在計算區(qū)域底部，初始計算時，使泥沙顆粒保持固定只計算由體積力驅(qū)動的水體。待計算流場紊流充分發(fā)展并穩(wěn)定后，使底床顆粒開始運動，再進行10 s的模擬。流場的初始剖面由對數(shù)律給定。為保證結(jié)果精度，取底床泥沙釋放后5～10 s的數(shù)據(jù)進行分析。

對于流體，x軸與y軸方向均設為循環(huán)邊界，頂邊界設為滑移邊界，底邊界為無滑移邊界。對于顆粒計算，計算區(qū)域四周均為循環(huán)邊界，頂邊界與底邊界為墻壁。

3 模擬結(jié)果分析

3.1 流體運動分析

此次模擬分為兩個階段。第一階段，水體由符合對數(shù)流速分布的初始剖面開始，受到體積力驅(qū)動，在固定泥沙顆粒的粗糙床面上運動，直至紊流發(fā)展充分。為觀察模擬流場的紊動特性，圖2和圖3給出了目標平均流速的算例在紊流充分發(fā)展后的流速剖面圖。圖2-a為接近泥沙顆粒頂部z=2.0 mm處的紊流流速剖面，圖2-b為z=2.0 mm與y=0.06 m兩平面交線處瞬時流速u和w沿y軸方向的分布。圖3-a、3-b兩圖分別為x=0 m與y=0.06 m處紊流流速的分布圖。

在接近泥沙顆粒頂部的z=2.0 mm軸剖面處，由圖2-a可以清楚地看到瞬時流速形成了低流速帶與高流速帶相間排列的結(jié)構(gòu)特征，與Kline等利用氫泡技術在光滑床面上觀察到的條帶結(jié)構(gòu)相似［14］。圖2-b中底床流速分量u和w沿y軸大小間隔分布的間距雖不盡相同，但其平均值卻相對穩(wěn)定，其無量綱化后距離約為424。

圖2 充分發(fā)展的近底紊流流速分布Fig.2 Distribution of fully developed turbulent flow profile near bottom

圖3 充分發(fā)展的側(cè)邊界紊流流速分布Fig.3 Distribution of fully developed turbulent flow profile at lateral boundary

從圖3-a中的x=0 m剖面可以看到，底部附近的低速流體被舉升抬離底床，發(fā)生“噴射”現(xiàn)象。圖3-b中進入主流區(qū)的低速流體受到高速流體的沖擊，產(chǎn)生大量紊動，低速流體被攜帶至下游方向，發(fā)生“清掃”過程。模擬中包含低速流體噴射與高速流體清掃的紊流猝發(fā)現(xiàn)象與學者對于此現(xiàn)象的一般認識吻合［15-16］。

模擬第二階段，固定在底床的顆粒由靜止狀態(tài)釋放，底床泥沙顆粒在受到水流作用后起動，隨著水流向下游運動。由于床面泥沙顆粒厚度不同，本文將9組算例按照初始泥沙床面高度分為小流速組0.2～0.6 m/s）和與大流速組。圖4中4-a與4-b為分別為和1.0 m/s兩組代表算例在第二階段的時均流速剖面圖，圖中的z坐標以計算域底邊界為零點。兩組代表算例在模擬區(qū)域底邊界上的流速均為零，說明底邊界上的泥沙顆粒并未運動，保證了泥沙厚度設置的合理性。圖4-a中的流速較小，在初始泥沙床面頂部高度以下的模擬流速值均接近零；而對于圖4-b中算例，初始床面頂部高度以下小范圍內(nèi)流速不為零。從以下3.2節(jié)泥沙運動的分析中可以得到，小流速組的算例泥沙運動形式為未起動或推移質(zhì)運動，泥沙均在初始床面附近運動，此時初始床面高度附近被泥沙占據(jù)，此處的流速接近于零；大流速組的泥沙運動為推移質(zhì)與懸移質(zhì)共存的形式，初始床面高度附近的泥沙被水流懸揚起來，因而此處的流速不為零，而且流速過渡較平滑。

初始床面以上的流速利用如下對數(shù)流速公式進行擬合

式中:κ為卡門常數(shù)取0.4，z1為理論床面距計算域底部高度，z為距計算域底部高度，ks為等效粗糙高度。對于m/s的算例，摩阻流速u*=0.045 m/s，理論床面高度z1=1.68 mm，等效粗糙高度ks=1.14D；對于的算例，摩阻流速u*=0.090 m/s，理論床面高度z1=4.80 mm，等效粗糙高度ks=3.13D。擬合結(jié)果如圖4-a與4-b所示，左上角小圖中z軸以理論床面為0點繪制，從圖可以看出平均流速均勻分布在理論曲線兩側(cè)，水流模擬結(jié)果是較為合理的。

圖4 時均流速剖面Fig.4 Time?averaged velocity profile

圖5 紊流強度分量與雷諾應力沿z軸方向的時間平均分布Fig.5 Distribution of relative turbulence intensity components and relative Reynold stress along the z axis

流體紊動強度的分布對于泥沙的運動有著至關重要的影響。圖5中的5-a、5-b、5-c圖為小流速代表組與大流速代表組按時間平均后的相對紊動強度分量沿z軸方向的分布圖。兩組紊動強度的峰值皆位于初始床面頂部附近，此處水流與泥沙相互作用最為劇烈。水流的紊動作用帶動泥沙的運動，泥沙的運動改變了床面的形狀，使得水流在床面附近的流動更為復雜，更容易產(chǎn)生紊動。在靠近模擬區(qū)域底邊界被泥沙大量占據(jù)，水流在孔隙中的流動受到較大限制，紊動強度接近于零。對于接近模擬區(qū)域頂部，流體沿x軸與y軸方向的流動受床面影響有限，故5-a、5-b中的紊流強度較小。流體沿z軸方向的速度在模擬區(qū)域頂部為零，故圖5-c中的紊動強度的軸分量為零。因為流體的流動方向為x軸正方向，紊動強度x軸的分量峰值較大，而y軸與z軸分量峰值接近。

圖5-d中為相對平均雷諾應力τR=-＜u'w'＞/u2

*沿z軸方向的分布。在計算區(qū)域上部，流體所受到的粘性應力相比于雷諾應力可以忽略，故雷諾應力應成線性分布，本文的模擬結(jié)果基本符合。在接近床面頂部處，泥沙從紊動中獲得能量從而運動，模擬的雷諾應力形成拐點。在模擬區(qū)域底邊界處，雷諾應力降為零。

圖6 模擬中三種典型的泥沙運動方式Fig.6 Three typical style of sediment transport in the simulations

圖7 泥沙濃度沿z軸方向的分布Fig.7 Profile of sediment concentration along z axis

3.2 泥沙運動分析

泥沙在不同的水流強度作用下有著不同的運動規(guī)律。選取目標流速＜＜u→f＞＞分別為0.2 m/s、0.6 m/s和1.0 m/s典型代表算例，將其流速分布與泥沙顆粒運動形式繪于圖6的6-a、6-b和6-c中。從中觀察得，目標流速為0.2 m/s的流場為層流狀態(tài)，處于床面的泥沙均處于靜止狀態(tài)，此時的泥沙顆粒并未起動；目標流速為0.6 m/s的流場底部有明顯的紊流結(jié)構(gòu)，床面泥沙顆粒隨水流沿著床面向下游運動，此時的泥沙為推移質(zhì)運動；目標流速為1.0 m/s的流場紊流結(jié)構(gòu)擴展至整個計算域，同時底床的泥沙有大量的起懸，泥沙為推移質(zhì)與懸移質(zhì)共同輸沙的狀態(tài)。由Shields曲線計算得到［1］，對于粒徑為0.5 mm的泥沙顆粒，其起動的臨界摩阻流速u*c為0.016 m/s，大于目標流速為0.2 m/s計算的摩阻流速0.012 m/s，驗證了泥沙運動模擬的合理性。根據(jù)泥沙的運動形態(tài)，本文的九組算例可以分為：未起動（目標流速0.2 m/s）、推移質(zhì)輸沙（目標流速0.3～0.6 m/s）和推移質(zhì)與懸移質(zhì)共同輸沙（目標流速0.7～1.0 m/s）三組。

以上三組的分組依據(jù)可以從計算區(qū)域內(nèi)泥沙濃度沿z軸方向的分布中得到。圖7繪制了9組算例的泥沙濃度分布圖，泥沙濃度沿z軸的分布范圍隨著目標流速的增大而增大。圖中目標流速為0.2 m/s的泥沙均集中于床面；目標流速為0.3～0.6 m/s的泥沙濃度隨著z軸高度的增加而減少，此組算例的泥沙床面無量綱高度為0.042 5，大量泥沙顆粒運動的范圍僅限于床面附近；目標流速為0.7～1.0 m/s的泥沙濃度分布則遍布整個z向計算范圍。在接近計算域底邊界處，泥沙的孔隙率為0.47，顆粒密為實堆積。對于推移質(zhì)與懸移質(zhì)共同輸沙組，在計算域頂邊界附近泥沙濃度發(fā)生振蕩。產(chǎn)生此種濃度振蕩原因為，泥沙與頂部壁面發(fā)生碰撞，運動被干擾，加之泥沙濃度較小，泥沙空間分布具有隨機性。

對于推移質(zhì)單寬輸沙率，前人已進行過大量的實驗并提出幾種較為合理的計算公式。為了驗證本文模擬的推移質(zhì)輸沙率強度，選取Gilbert、Meyer?Peter、Wilson及錢寧的實驗值，Einstein、Meyer?Peter、Yalin及En?gelund經(jīng)驗計算公式［17］，與本文目標流速為0.3～0.6m/s算例的統(tǒng)計值進行比較。對比結(jié)果如圖8所示，橫軸為推移質(zhì)輸沙強度Φ，縱軸為Shields數(shù)Ψ（水流強度的倒數(shù)），兩者的計算公式如下

式中：gb為推移質(zhì)單寬輸沙率；D為泥沙粒徑；γp和γf分別為顆粒和流體的重度；τ0為作用在單位床面面積上的水流拖曳力。

從圖8可以看到，本文推移質(zhì)輸沙率的模擬值與Gilbert、Meyer?Peter等實驗值及四種經(jīng)驗公式比較接近，這驗證了模型對于水流作用下顆粒運動描述的合理性。

圖8 推移質(zhì)輸沙率的模擬值、實驗值與理論公式比較Fig.8 Comparison of bedload transportation among simulation result,experimental data and analytical formula

4 結(jié)論

本文采用CFD?DEM方法，研究了九組不同水動力條件下流體與泥沙顆粒的運動狀態(tài)。模擬中，流體計算采用Smagorinsky紊流模型，顆粒接觸計算采用Hertz?Mindlin無滑移“軟球”模型，利用開源軟件CFDEM耦合庫進行流體-顆粒作用力的交換，完成耦合過程。從水流與顆粒運動兩方面分析模擬結(jié)果得到以下結(jié)論：

（1）水流在床面附近形成明顯的高低流速間隔的條帶狀結(jié)構(gòu)，并且可以觀察到紊流猝發(fā)現(xiàn)象，符合現(xiàn)有對紊流結(jié)構(gòu)的認識。

（2）模擬得到的流速平均剖面符合對數(shù)律公式，保證了水動力條件的合理性。

（3）分析得到紊動強度三個方向分量與雷諾應力沿z軸的分布，紊動強度峰值出現(xiàn)在床面附近，接近底床和計算頂邊界強度較小。計算區(qū)域上部流體的雷諾應力與理論分布接近，成線性變化。

（4）9組泥沙顆粒運動模擬中包括未起動、推移質(zhì)輸沙及推移質(zhì)與懸移質(zhì)共同輸沙，平均泥沙濃度垂向分布表明了三種狀態(tài)的存在。

（5）模擬得到的推移質(zhì)輸沙率與實驗值及理論公式吻合良好。

［1］錢寧，萬兆恵.泥沙運動力學［M］.北京：科學出版社，1983.

［2］王光謙，胡春宏.泥沙研究進展［M］.北京：中國水利水電出版社，2006.

［3］王光謙，孫其誠.顆粒物質(zhì)及其多尺度結(jié)構(gòu)統(tǒng)計規(guī)律［J］.工程力學，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.

［4］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.

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

［6］Ji C N，Ante M，Eldad A，et al.Numerical investigation of particle saltation in the bed?load regime［J］.Science China Technologi?cal Sciences，2014，57（8）：1 500-1 511.

［7］及春寧，劉丹青，許棟.基于顆粒離散元的沙紋演變大渦模擬研究［J］.力學學報，2015，47（4）：613-623. JI C N，LIU D Q，XU D.Large eddy simulation of sand ripple evolution using Discrete Particle Method［J］.Chinese Journal of Theo?retical and Applied Mechanics，2015，47（4）：613-623.

［8］Pan T W，Joseph D D，Bai R，et al.Fluidization of 1204 spheres：simulation and experiment［J］.Journal of Fluid Mechanics，2002，451：169-191.

［9］Ji C，Munjiza A，Williams J J R.A novel iterative direct?forcing immersed boundary method and its finite volume applications［J］. Journal of Computational Physics，2012，231（4）：1 797-1 821.

［10］喬光全，張慶河，張金鳳.扁球體泥沙顆粒沉降的三維格子Boltzmann模型［J］.泥沙研究，2013（4）：1-7. QIAO G Q，ZHANG Q H，ZHANG J F.Three dimensional Lattice Boltzmann model for settling of oblate spheroid sediment parti?cles［J］.Journal of Sediment Research，2013（4）：1-7.

［11］Kloss C，Goniva C，Hager A，et al.Models，algorithms and validation for opensource DEM and CFD?DEM［J］.Progress in Computa?tional Fluid Dynamics，an International Journal，2012，12（2-3）：140-152.

［12］Smagorinsky J.General circulation experiments with the primitive equations：I.the basic experiment［J］.Monthly weather review，1963，91（3）：99-164.

［13］孫其誠，王光謙.顆粒物質(zhì)力學導論［M］.北京：科學出版社，2009.

［14］Kline S J，Reynolds W C，Schraub F A，et al.The structure of turbulent boundary layers［J］.Journal of Fluid Mechanics，1967，30 （4）：741-773.

［15］Grass A J.Structural features of turbulent flow over smooth and rough boundaries［J］.Journal of Fluid Mechanics，1971，50（2）：233-255.

［16］建忠.湍流的擬序結(jié)構(gòu)［M］.北京：機械工業(yè)出版社，1995.

［17］錢寧.推移質(zhì)公式的比較［J］.水利學報，1980（4）：1-11. QIAN N.A Comparison of the Bed Load Formulas［J］.Journal of Hydraulic Engineering，1980（4）：1-11.

Simulation of fluid?sediment particle motion near bed based on CFD?DEM coupling method

SU Dong?sheng,ZHANG Qing?he,SUN Jian?jun,LI Ming?xing
（State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin 300072,China）

Based on CFD?DEM method,using computational fluid dynamics software OpenFOAM,particle mo?tion simulation software LIGGGHTS and CFDEM coupling library,a coupled fluid?particle model was developed and applied in the investigation of flow and particle near bed in different hydrodynamic conditions.The simulation results show that the coupled model has a good performance on describing turbulence intensity and the distribution of Reynold stress.For particle motion,this model can also depict well the three typical states of settled,bedload transport and bedload?suspended load transport.The simulated results of bedload discharge per unit have a good match with experimental data and empirical formula.

CFD?DEM coupled model;near bed flow;turbulence intensity;bedload sediment transport

TV 142；O 242.1

A

1005-8443（2016）03-0224-07

2015-11-16；

2015-11-25

國家自然科學基金資助（51179122）

蘇東升（1991-），男，陜西省寶雞人，碩士研究生，主要從事港口海岸及近海工程研究。

Biography：SU Dong?sheng（1991-），male，master student.