光滑壁面明渠水流的三維格子玻耳茲曼模擬

丁 磊，張慶河

(天津大學(xué)建筑工程學(xué)院港口與海洋工程教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072)

明渠流是自然界和工程實(shí)踐中常見的一種水流運(yùn)動，泥沙運(yùn)動以及河床演變等問題都與之密切相關(guān)．國內(nèi)外學(xué)者主要采用實(shí)驗(yàn)、理論分析和數(shù)值模擬的方法進(jìn)行明渠均勻流的研究．在實(shí)驗(yàn)研究方面，Nezu等[1]采用二維激光多普勒流速儀(LDA)研究了水槽中光滑壁面均勻紊流，認(rèn)為近底附近的時均流速分布可以用對數(shù)律來描述，同時還存在著分區(qū)結(jié)構(gòu)，后續(xù)的相關(guān)實(shí)驗(yàn)[2-3]也證實(shí)了這一觀點(diǎn)．實(shí)驗(yàn)研究為認(rèn)識明渠均勻流的水流特性和紊動結(jié)構(gòu)奠定了基礎(chǔ)．在理論分析方面，研究者主要從恒定明渠紊流的雷諾方程出發(fā)推導(dǎo)光滑壁面時均流速、切應(yīng)力以及雷諾應(yīng)力分布[4-6]．在明渠紊流的數(shù)值模擬方面，各種紊流模型、大渦模擬(LES)以及直接數(shù)值模擬(DNS)都得到了應(yīng)用[7-9]．?dāng)?shù)值模擬結(jié)果提供了比實(shí)驗(yàn)更為詳細(xì)的明渠水流紊動特性和整體流場信息，為進(jìn)一步研究明渠流中的泥沙運(yùn)動規(guī)律打下了基礎(chǔ)．但由于問題的復(fù)雜性，在水流與泥沙顆粒耦合模擬方面還有很多問題有待深入研究．為此，尋找新型數(shù)值模擬方法研究明渠水流，特別是明渠水流作用下的泥沙運(yùn)動仍是值得探索的工作．

近年來，格子玻耳茲曼(lattice Boltzmann，LB)方法在許多學(xué)科領(lǐng)域得到了廣泛應(yīng)用．LB方法具有易于處理固體邊界、易于并行計(jì)算等特點(diǎn)，較傳統(tǒng)的數(shù)值計(jì)算方法更容易處理動邊界且更容易對流體動力作用下的大規(guī)模顆粒運(yùn)動進(jìn)行全尺度直接模擬[10-13]．因此，LB方法研究明渠流中的泥沙運(yùn)動具有較大的優(yōu)勢．目前，LB方法在明渠水流中的應(yīng)用主要集中于二維模型[14-16]，張金鳳等[17]、Fernandino 等[18]采用三維LB結(jié)合大渦模擬(LES)對光滑明渠流進(jìn)行了模擬，但所得到的紊流時均流速分布與解析解或?qū)嶒?yàn)值還有一定差距．由于 LES對紊流的模擬存在一定程度的簡化，特別是文獻(xiàn)中采用的 LES主要是傳統(tǒng)的Smagorinsky模式[17-18]，而相關(guān)研究表明該模式對壁面附近紊流的描述存在缺陷[19]．鑒于此，筆者進(jìn)一步采用三維 LB方法對光滑壁面明渠紊流進(jìn)行直接數(shù)值模擬(DNS)，并與解析解和實(shí)驗(yàn)結(jié)果進(jìn)行比較，以便為應(yīng)用 LB方法研究明渠流中泥沙運(yùn)動機(jī)理奠定基礎(chǔ)．

1 數(shù)值方法

1.1 三維LB方法

與傳統(tǒng)的計(jì)算流體力學(xué)方法直接離散求解Navier-Stokes(N-S)方程不同，LB方法是在介觀層面上構(gòu)建簡化的離散速度動力學(xué)模型，即格子模型，通過控制格子模型中虛擬粒子的密度分布函數(shù)演化過程使計(jì)算域中流體在宏觀層面上表現(xiàn)出 N-S方程控制的性質(zhì)．演化過程的控制方程即為 LB方程[20]，如式(1)所示．

式中：r為位置矢量；ei為格子模型中離散速度矢量；i為離散速度方向；t為時間變量；tΔ為時間步長；if為密度分布函數(shù)；iΔ為碰撞算子．

密度分布函數(shù)與宏觀流體變量(密度ρ和動量密度ρu)的關(guān)系如式(2)和式(3)所示．

在本文的三維計(jì)算中，選取了D3Q19格子模型，即離散速度的個數(shù)為 19，離散速度矢量的具體分布形式如式(4)所示．

1.2 LB模擬流動的驅(qū)動方式

LB中常用的產(chǎn)生流體運(yùn)動的驅(qū)動方式主要有3種[21]：①將流體所受體積力作為外力項(xiàng)加入 LB方程來實(shí)現(xiàn)；②作用壓力梯度，而對于不可壓縮流體則往往將壓力梯度轉(zhuǎn)化成體積力作為驅(qū)動力；③在入口邊界給定速度分布．張金鳳等[17]采用體積力的方式實(shí)現(xiàn)明渠水流的驅(qū)動，但給出的體積力的計(jì)算方法主要適用于層流情況，在明渠紊流中使用則會由于雷諾應(yīng)力的影響帶來一定的偏差．Fernandino等[18]也將壓力梯度轉(zhuǎn)化成體積力來進(jìn)行計(jì)算，所得到的時均流速分布與實(shí)驗(yàn)值同樣存在差距．為消除壓力作為體積力驅(qū)動時可能形成的誤差，這里筆者采用入口速度驅(qū)動的方式．

1.3 邊界條件

圖 1給出了模擬時的三維坐標(biāo)、水流流動方向以及各邊界的位置．計(jì)算區(qū)域的入口設(shè)置速度邊界條件[22]，側(cè)邊界和出口處施加周期性邊界條件，底面采用無滑移邊界條件，自由表面則采用自由滑移邊界條件[21]進(jìn)行近似處理．由于這里主要研究明渠近底水動力特性，選取自由滑移邊界條件是合適的．

圖1 計(jì)算區(qū)域與邊界示意Fig.1 Sketch of computational domain and boundary

1.4 紊流的描述

明渠底面的存在使得水流紊動變得更加復(fù)雜．底面附近紊動的合理描述對認(rèn)識明渠流近底水動力特性以及揭示泥沙顆粒運(yùn)動機(jī)理具有重要意義．在數(shù)值模擬中，由于DNS無需對流動做任何簡化或近似，可以比較完整準(zhǔn)確地反映紊動特性，已經(jīng)越來越多地應(yīng)用于紊流特別是壁面紊流的研究中[23-24]．因此，本文采用DNS方法來模擬近底水流紊動．

2 明渠均勻流模擬

2.1 模型驗(yàn)證

為驗(yàn)證建立的模型，首先計(jì)算明渠層流情況．明渠層流流速沿水深的分布如式(5)所示，滿足拋物線關(guān)系[25]．

式中：u( y)為距底面距離為y處的流速值；umax為流速分布的最大值；h為明渠水深．這里取h＝1.0×10-2,m，umax＝7.2×10-2,m/s，以斷面平均流速 Um定義雷諾數(shù)Re＝Umh/ν＝480，v為水的運(yùn)動黏滯系數(shù)．計(jì)算區(qū)域沿流向、水深及寬度方向分別取為10×10-2、1×10-2和2×10-2,m．需要說明的是，由于流動處于層流狀態(tài)時，各層流體之間相互影響較小，因此為減少計(jì)算時間，在水深方向上只選取了底面以上1×10-2,m的區(qū)域．采用立方體網(wǎng)格離散計(jì)算域，沿流向、水深及寬度方向網(wǎng)格數(shù)分別為1,000、100、200，網(wǎng)格大小為Δx=Δy=Δz = 1× 1 0-4m．計(jì)算中將式(5)作為入口速度邊界條件．

明渠層流中各層切應(yīng)力τ,如式(6)所示，容易看出，切應(yīng)力沿水深呈線性分布．

式中μ為水的動力黏滯系數(shù)．

圖2給出了計(jì)算區(qū)域中間斷面流速分布與式(5)的比較結(jié)果，圖3給出了切應(yīng)力分布與式(6)的比較結(jié)果．由圖2和圖3可知，流速分布與切應(yīng)力分布的模擬結(jié)果與解析解有很好的一致性，這說明所建立的 LB模型可以應(yīng)用于低雷諾數(shù)時明渠均勻?qū)恿鞯哪M．

圖2 明渠層流流速分布Fig.2 Velocity profile of laminar open channel flow

圖3 明渠層流切應(yīng)力分布Fig.3 Shear stress profile of laminar open channel flow

2.2 光滑壁面明渠紊流近底水流模擬

2.2.1 算例設(shè)置

進(jìn)一步將 LB模型應(yīng)用于光滑壁面明渠紊流近底運(yùn)動的模擬，為便于與實(shí)驗(yàn)結(jié)果比較，計(jì)算時選取Nezu等[1]的實(shí)驗(yàn)條件，具體參數(shù)見表1．

表1 實(shí)驗(yàn)參數(shù)Tab.1 Experimental parameters

根據(jù)以往實(shí)驗(yàn)研究和理論分析的成果[1-2]，光滑壁面明渠紊流近底存在著時均流速分布的分區(qū)結(jié)構(gòu)，沿水深方向由下至上依次為黏性底層、過渡層和對數(shù)層．黏性底層內(nèi)的流速呈線性分布，對數(shù)層內(nèi)的流速符合對數(shù)律，過渡層由于受邊界條件等因素影響較大，尚沒有統(tǒng)一的流速分布表達(dá)式．在設(shè)置入口速度邊界條件時考慮了分區(qū)結(jié)構(gòu)，其中黏性底層和對數(shù)層中的流速分布表達(dá)式均引自文獻(xiàn)[1]，這里對過渡層的流速分布做了簡化處理，在已知黏性底層和對數(shù)層流速后利用二次曲線擬合連接，具體形式見式(7)．

式中：u+=u/ u*1；y+= u*1y/ν．為了加快紊動的發(fā)展，入口邊界條件在式(7)的基礎(chǔ)上疊加了隨機(jī)小幅擾動．計(jì)算區(qū)域選擇在實(shí)驗(yàn)渠道的中心斷面附近，沿流向、水深方向、寬度方向分別取10×10-2,m、3×10-2,m、2×10-2,m．需要指出的是，由于主要研究近底水動力條件，考慮到DNS的計(jì)算較為耗時，在保證計(jì)算精度和計(jì)算效率的情況下，本文在水深方向減小了計(jì)算區(qū)域，只選取了底面以上3×10-2,m 的區(qū)域，此時計(jì)算域的上邊界仍處于對數(shù)層，受自由表面影響很?。?jì)算域采用立方體網(wǎng)格離散，沿流向、水深和寬度方向的網(wǎng)格數(shù)分別為 1,000、300和 200，網(wǎng)格大小為Δx = Δy = Δz = 1× 1 0-4m，換算至量綱1的壁面單位表示為 Δ x+=Δ y+=Δ z+= u*1Δ x /ν= 0 .431，網(wǎng)格精度滿足DNS的計(jì)算要求[26]．

2.2.2 流速分布

圖4(a)顯示了計(jì)算區(qū)域中部第60,000時間步的瞬時流速垂向剖面，圖 4(b)則給出了時均流速垂向剖面(取最后 40,000時間步進(jìn)行時間平均)與式(7)和實(shí)驗(yàn)結(jié)果的比較．從瞬時流速分布看，流速剖面較好地反映了垂向分區(qū)特性．由圖 4(b)可知，計(jì)算值與理論分析及實(shí)驗(yàn)結(jié)果均符合較好，說明式(7)中采用二次曲線擬合過渡層的流速分布是合理的．

圖4 光滑壁面明渠紊流流速分布Fig.4 Velocity distribution of smooth turbulent open channel flow

2.2.3 切應(yīng)力

總的時均切應(yīng)力由時均黏性應(yīng)力和雷諾應(yīng)力兩部分構(gòu)成．圖5給出了明渠中部沿流向剖面的無量綱時均黏性應(yīng)力τ1*= (μ? u / ? y ) /的等值線分布．由圖5可知，時均黏性應(yīng)力在壁面附近最大，在距壁面一定位置后迅速減小，這與已有的認(rèn)識是一致的[1,27]．

圖5 時均黏性應(yīng)力等值線分布Fig.5 Contours of time-averaged viscous stress

無量綱雷諾應(yīng)力可用τ*=-'/ u2表示，圖 6給出了τ2*沿水深分布的模擬結(jié)果與實(shí)驗(yàn)結(jié)果的比較情況．兩者最大相對誤差的絕對值為 6.8%，說明計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合較好．同時圖 6還反映了雷諾應(yīng)力沿水深方向除壁面附近外基本為線性分布，這也與實(shí)驗(yàn)和理論分析結(jié)果相符[1,27]．

底面時均切應(yīng)力是流動阻力規(guī)律研究的關(guān)鍵，特別是由其定義的摩阻流速對流速分布和紊動特性有重要影響．利用計(jì)算結(jié)果并根據(jù)董曾南等[2]給出的方法計(jì)算了底面時均切應(yīng)力，并將其換算成摩阻流速(0.389×10-2,m/s)，與實(shí)驗(yàn)得到的摩阻流速(0.431×10-2m/s)相比，相對誤差為 9.74%，說明計(jì)算結(jié)果可以較好地反映底面時均切應(yīng)力．由于計(jì)算所采用的實(shí)驗(yàn)條件中明渠水流的寬深比為 5.9，中心區(qū)域受側(cè)壁影響較小，底面時均切應(yīng)力沿寬度方向即橫向分布應(yīng)較為均勻，圖7給出的中間斷面底面時均切應(yīng)力τ0的橫向分布說明了這一點(diǎn)．

圖6 雷諾應(yīng)力分布Fig.6 Reynolds stress distribution

圖7 底面時均切應(yīng)力橫向分布模擬結(jié)果Fig.7 Lateral distribution of bottom shear stress

3 結(jié) 語

應(yīng)用三維 LB方法對明渠均勻流進(jìn)行了模擬，為了避免以往 LB計(jì)算中壓力作為體積力驅(qū)動水流可能形成的誤差，采用了入口速度驅(qū)動方式．將該方式應(yīng)用于明渠層流模擬時，流速和切應(yīng)力沿水深分布與解析解吻合很好．對于光滑壁面明渠紊流，采用 DNS方式模擬水體紊動，并重點(diǎn)對近底水動力特性進(jìn)行了模擬研究．模擬得到的時均流速分布與考慮水流垂向分區(qū)結(jié)構(gòu)的流速分布理論解和實(shí)驗(yàn)結(jié)果均符合較好．另外，雷諾應(yīng)力的垂向剖面與實(shí)驗(yàn)結(jié)果有較好的一致性，從底面時均切應(yīng)力的橫向分布結(jié)果也可以看出，建立的模型可以較好地描述近底紊動特性和阻力規(guī)律．

從有限的算例結(jié)果看，采用的設(shè)置入口速度邊界條件并結(jié)合直接數(shù)值模擬的 LB方法，可較準(zhǔn)確地反映光滑明渠均勻?qū)恿髋c紊流的近底水動力特性．在此基礎(chǔ)上，研究不同雷諾數(shù)下光滑明渠近底紊流的水動力特性，并進(jìn)一步研究粗糙壁面明渠水流的近底水動力特性，從微觀角度探索泥沙顆粒運(yùn)動機(jī)理，是下一步需要深入的工作．

［1］Nezu I，Rodi W. Open-channel flow measurements with a laser Doppler anemometer[J].Journal of Hydraulic Engineering，1986，112(5)：335-355.

［2］董曾南，丁 元. 光滑壁面明渠均勻紊流水力特性[J].中國科學(xué)A輯，1989(11)：1208-1218.

Dong Zengnan，Ding Yuan. Turbulence characteristics of uniform flow in a smooth open channel[J].Science in China，Ser A，1989(11)：1208-1218(in Chinese).

［3］劉春晶，李丹勛，王興奎. 明渠均勻流的摩阻流速及流速分布[J]. 水利學(xué)報(bào)，2005，36(8)：950-955.

Liu Chunjing，Li Danxun，Wang Xingkui. Experimental study on friction velocity and velocity profile of open channel flow[J].Journal of Hydraulic Engineering，2005，36(8)：950-955(in Chinese).

［4］Yang S-Q，McCorquodale J A. Determination of boundary shear stress and Reynolds shear stress in smooth rectangular channel flows[J].Journal of Hydraulic Engineering，2004，130(5)：458-462.

［5］Guo J，Julien P Y. Shear stress in smooth rectangular open-channel flows[J].Journal of Hydraulic Engineering，2005，131(1)：30-37.

［6］劉亞坤，倪漢根. 水力光滑明渠流流速分布的新公式[J]. 水利學(xué)報(bào)，2007，38(11)：1336-1340.

Liu Yakun，Ni Hangen. New formula for velocity profile of turbulent flow in open channel with smooth wall[J].Journal of Hydraulic Engineering，2007，38(11)：1336-1340(in Chinese).

［7］Kang H，Choi S-U. Reynolds stress modelling of rectangular open-channel flow[J].International Journal for Numerical Methods in Fluids，2006，51(11)：1319-1334.

［8］Hinterberger C，F(xiàn)r?hlich J，Rodi W. 2D and 3D turbulent fluctuations in open channel flow withReτ＝590 studied by large eddy simulation[J].Flow，Turbulence and Combustion，2008，80(2)：225-253.

［9］Honda I，Kanagawa N，Kawashima Y. Direct numerical simulation of open channel flow with buoyancy[J].Developments in Chemical Engineering and Mineral Processing，2003，11(5/6)：465-479.

［10］Ladd A J C，Verberg R. Lattice-Boltzmann simulations of particle-fluid suspensions[J].Journal of Statistical Physics，2001，104(5/6)：1191-1251.

［11］吳錘結(jié)，周菊光. 懸浮顆粒運(yùn)動的格子Boltzmann數(shù)值模擬[J]. 力學(xué)學(xué)報(bào)，2004，36(2)：151-162.

Wu Chuijie，Zhou Juguang. Numerical simulation of suspension motion of irregular shaped particles via the lattice Boltzmann method[J].Acta Mechanica Sinica，2004，36(2)：151-162(in Chinese).

［12］Stratford K，Pagonabarraga I. Parallel simulation of particle suspensions with the lattice Boltzmann method[J].Computers and Mathematics with Applications，2008，55(7)：1585-1593.

［13］H?lzer A，Sommerfeld M. Lattice Boltzmann simulations to determine drag，lift and torque acting on non-spherical particles[J].Computers and Fluids，2009，38(3)：572-589.

［14］程永光，索麗生. 二維明渠非恒定流的格子 Boltzmann模擬[J]. 水科學(xué)進(jìn)展，2003，14(1)：9-14.

Cheng Yongguang，Suo Lisheng. 2-D open channel flow simulations by the lattice Boltzmann model[J].Advances in Water Science，2003，14(1)：9-14(in Chinese).

［15］Li Y，Huang P. A coupled lattice Boltzmann model for the shallow water-contamination system[J].International Journal for Numerical Methods in Fluids，2008，59(2)：195-213.

［16］Liu H，Zhou J G，Burrows R. Numerical modeling of turbulent compound channel flow using the lattice Boltzmann method[J].International Journal for Numerical Methods in Fluids，2009，59(7)：753-765.

［17］張金鳳，張慶河. 明渠流中分形絮團(tuán)破裂的格子Boltzmann方法模擬[J]. 天津大學(xué)學(xué)報(bào)，2009，42(1)：17-23.

Zhang Jinfeng，Zhang Qinghe. Lattice Boltzmann simulation of the breakup process of fractal flocs in openchannel flows[J].Journal of Tianjin University，2009，42(1)：17-23(in Chinese).

［18］Fernandino M，Beronov K，Ytrehus T. Large eddy simulation of turbulent open duct flow using a lattice Boltzmann approach[J].Mathematics and Computers in Simulation，2009，79(5)：1520-1526.

［19］Sagaut P.Large Eddy Simulation for Incompressible Flows：An Introduction[M]. 3rd ed. Berlin：Springer，2006.

［20］Chen S，Doolen G D. Lattice Boltzmann method for fluid flows[J].Annual Review of Fluid Mechanics，1998，30(1)：329-364.

［21］Dupuis A. From a Lattice Boltzmann Model to a Parallel and Reusable Implementation of a Virtual River[D]. Geneva：University of Geneva，2002.

［22］He X Y，Zou Q S，Luo L S，et al. Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model[J].Journal of Statistical Physics，1997，87(1/2)：115-136.

［23］Lammers P，Beronov K N，Volkert R，et al. Lattice BGK direct numerical simulation of fully developed turbulence in incompressible plane channel flow[J].Computers and Fluids，2006，35(10)：1137-1153.

［24］Wu X，Moin P. Direct numerical simulation of turbulence in a nominally zero-pressure-gradient flat-plate boundary layer[J].Journal of Fluid Mechanics，2009，630：5-41.

［25］格拉夫，阿廷拉卡. 河川水力學(xué)[M]. 趙文謙，萬兆惠，譯. 成都：成都科技大學(xué)出版社，1997.

Graf W H，Altinakar M S.River Hydraulics[M]. Zhao Wenqian，Wan Zhaohui，Trans. Chengdu：Chengdu Science and Technology University Press，1997(in Chinese).

［26］Singh K M，Sandham N D，Williams J J R. Numerical simulation of flow over a rough bed[J].Journal of Hydraulic Engineering，2007，133(4)：386-398.

［27］章梓雄，董曾南. 粘性流體力學(xué)[M]. 北京：清華大學(xué)出版社，1998.

Chwang A T，Dong Zengnan.Viscous Fluid Mechanics[M]. Beijing：Tsinghua University Press，1998(in Chinese).