沙紋床面振蕩流邊界層的三維格子玻爾茲曼模擬研究

周志博 ，張慶河，張金鳳

(1.天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2.交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究院海洋水動(dòng)力研究中心，天津 300456)

沙質(zhì)海岸海床床面在波浪作用下經(jīng)常呈現(xiàn)沙紋形態(tài)，沙紋的出現(xiàn)會(huì)顯著增加床面對(duì)底部水體的阻力，使波浪底部邊界層流動(dòng)產(chǎn)生明顯改變．因此，對(duì)沙紋床面振蕩流邊界層流動(dòng)特性的深入認(rèn)識(shí)，是研究存在沙紋情況下的泥沙運(yùn)動(dòng)規(guī)律、準(zhǔn)確計(jì)算輸沙量等問(wèn)題的前提．為此，不少學(xué)者通過(guò)實(shí)驗(yàn)、數(shù)值模擬和現(xiàn)場(chǎng)觀測(cè)等手段對(duì)沙紋床面振蕩流邊界層的水流及其紊動(dòng)特性進(jìn)行了研究．實(shí)驗(yàn)研究方面，學(xué)者們通過(guò)激光流速儀(laser Doppler anemometry，LDA)[1]、多普勒剖面流速儀[2]等設(shè)備測(cè)量了沙紋床面邊界層流速、紊動(dòng)強(qiáng)度和床面切應(yīng)力等，為認(rèn)識(shí)沙紋床面振蕩流邊界層流動(dòng)奠定了基礎(chǔ)．在數(shù)值模擬方面，受到計(jì)算資源的制約，早期的工作主要進(jìn)行了二維數(shù)值模擬研究[1,3-7]．由于沙紋床面振蕩流邊界層內(nèi)的流動(dòng)通常為高雷諾數(shù)三維紊流，只有通過(guò)三維數(shù)值模擬才能較全面、細(xì)致地反映復(fù)雜的振蕩流邊界層流動(dòng)，為此，近年來(lái)人們開(kāi)始基于 DNS[8]、LES[9]等方法開(kāi)展三維數(shù)值模擬研究，以加強(qiáng)對(duì)沙紋邊界層復(fù)雜紊動(dòng)特性的認(rèn)識(shí)．然而，由于沙紋振蕩流邊界層流動(dòng)的復(fù)雜性，對(duì)沙紋存在導(dǎo)致的流動(dòng)分離、渦旋發(fā)展，沙紋床面切應(yīng)力的變化過(guò)程和沙紋形狀摩阻與膚面摩阻在總拖曳力中的貢獻(xiàn)等問(wèn)題仍有待深入研究．

格子玻爾茲曼(lattice Boltzmann，LB)方法是近年來(lái)發(fā)展迅速的一種新型數(shù)值模擬方法，可以從介觀角度描述宏觀流體的運(yùn)動(dòng)，具有可以高效并行計(jì)算、易于處理復(fù)雜幾何邊界等優(yōu)點(diǎn)．二維LB模型對(duì)光滑平床振蕩流邊界層獲得了較好的模擬結(jié)果[10-12]．對(duì)于粗糙床面振蕩流問(wèn)題；丁磊等[13]采用三維LB方法研究了振蕩流邊界層中顆粒的受力，并對(duì)球形顆粒規(guī)則排列構(gòu)成的粗糙平床上振蕩流邊界層流動(dòng)進(jìn)行了三維全尺度 DNS模擬，研究了粗糙床面的阻力特性[14]；周志博等[15]建立了基于Smagorinsky大渦模擬方法的三維LB模型，研究了粗糙床面振蕩流邊界層的紊動(dòng)特性和床面粗糙參數(shù)．

然而對(duì)于紊動(dòng)更加復(fù)雜的沙紋床面振蕩流邊界層流動(dòng)，通常振幅雷諾數(shù)較高，采用 DNS求解計(jì)算量巨大，而傳統(tǒng)的LES方法，如Smagorinsky模型則不能完全滿足壁面律[16]．為了能更好地反映近壁紊動(dòng)特性，學(xué)者們[16]提出了一種新的壁面自適應(yīng)局部渦黏(wall-adapting local eddy-viscosity，WALE)模型，該模型的建立不需要給定距離壁面的距離，也不需要引入衰減函數(shù)或求解動(dòng)力過(guò)程，并且能合理地描述轉(zhuǎn)捩流態(tài)的特性，更適用于復(fù)雜邊界情況[16-17]．因此本文將 WALE模型引入到三維 LB模型，以更好地描述沙紋振蕩流邊界層流動(dòng)．

1 數(shù)值方法

1.1 格子玻爾茲曼方法

格子玻爾茲曼是在介觀層面構(gòu)建離散速度模型，即格子模型，通過(guò)控制虛擬粒子在格子模型中的演化使流體表現(xiàn)出宏觀運(yùn)動(dòng)方程控制的性質(zhì)．LB方法的控制方程，即LB方程為

式中fi(r,t)為t時(shí)刻在r位置處沿i方向的密度分布函數(shù)；ei為沿i方向的離散速度；Di為碰撞算子，表示虛擬粒子的碰撞規(guī)則．

根據(jù)碰撞算子的選擇不同，可以將LB模型分為單松弛時(shí)間(single-relaxation-time，SRT)模型[18-19]和多松弛時(shí)間(multiple-relaxation-time，MRT)模型[20-21].LBGK模型[18-19]是 SRT模型中最具代表性的一種模型.Ladd等[22]在研究顆粒流問(wèn)題時(shí)發(fā)現(xiàn)，采用LBGK模型的計(jì)算結(jié)果在流固邊界處容易產(chǎn)生誤差，認(rèn)為這是在整個(gè)計(jì)算域采用同一個(gè)松弛時(shí)間造成的，并對(duì)LBGK模型進(jìn)行了改進(jìn)，提出了雙松弛時(shí)間(tworelaxation-time，TRT)模型．該模型的控制方程可以表示為

密度分布函數(shù)與流體宏觀物理量，即流體密度r、動(dòng)量密度j=ru和動(dòng)量通量Π之間的關(guān)系為

式中u為流體宏觀速度矢量．

離散速度模型是反映虛擬粒子速度的離散模式.本文采用D3Q19模型，各方向離散速度表達(dá)式為

1.2 邊界條件

1.2.1 周期性邊界條件

在LB模型中，周期性邊界上未知的分布函數(shù)可以通過(guò)虛擬粒子在格子模型中的運(yùn)動(dòng)規(guī)則直接確定，當(dāng)代表流體的虛擬粒子從一側(cè)周期性邊界離開(kāi)計(jì)算區(qū)域時(shí)，在下一個(gè)時(shí)間步會(huì)從計(jì)算區(qū)域另一側(cè)的周期性邊界進(jìn)入流場(chǎng)．以x方向采用周期性邊界為例，計(jì)算域在x方向劃分Nx個(gè)網(wǎng)格，該格式可表示為

式中fi為未知分布函數(shù)，其離散速度ei方向指向流場(chǎng)內(nèi)部．

1.2.2 自由滑移邊界條件

采用鏡面反射格式處理 LB模型中的自由滑移邊界．當(dāng)粒子由流體格點(diǎn)遷移至邊界格點(diǎn)后，沿著邊界法線方向的對(duì)稱方向反射．該格式可以表示為

式中：ei￠為ei關(guān)于邊界法向方向鏡像對(duì)稱的離散速度；rb和rf分別表示邊界格點(diǎn)和流體格點(diǎn)的位置，rb=rf+eiDt．鏡面反射格式的粒子在自由滑移邊界上不進(jìn)行碰撞，只進(jìn)行分布函數(shù)的遷移．按鏡面反射格式處理的邊界法向速度分量為0．

1.2.3 固壁邊界條件

采用具有二階精度的 Half-way反彈格式處理沙紋表面的固壁邊界[23]．Half-way反彈格式的邊界位于相鄰的固體格點(diǎn)與流體格點(diǎn)的中心處，即(rb+rf) /2處．具體的處理方式表示為

式中表示i的反方向．Half-way反彈格式的物理圖像為：在t時(shí)刻，流體節(jié)點(diǎn)處經(jīng)過(guò)碰撞后速度指向邊界的粒子進(jìn)行遷移，經(jīng)過(guò)Dt/2后到達(dá)處的邊界，然后執(zhí)行反彈，再經(jīng)過(guò)Dt/2后重新到達(dá)流體格點(diǎn)位置．圖 1為 Half-way反彈格式的曲線邊界示意[24]，為方便描述，這里以二維情況說(shuō)明．圖中方形格點(diǎn)為采用 Half-way反彈格式的邊界格點(diǎn)，位于流體格點(diǎn)和固體格點(diǎn)之間，黑色實(shí)線所示為實(shí)際邊界，模型中通過(guò)空心方塊連接而成的邊界近似表示.

圖1 Half-way反彈格式曲線邊界示意Fig.1 Curved boundary of Half-way bounce-back boundary condition

此外，固壁邊界受到流體的作用力可以通過(guò)動(dòng)量交換法[25]計(jì)算．對(duì)某固壁邊界點(diǎn)，相鄰的流體格點(diǎn)向邊界方向遷移的分布函數(shù)為，對(duì)應(yīng)的動(dòng)量為；經(jīng)過(guò)反彈后離開(kāi)固壁邊界點(diǎn)的分布函數(shù)為，對(duì)應(yīng)的動(dòng)量為；由動(dòng)量定理和力的相互作用原理可知，固壁邊界點(diǎn)受到的作用力大小為；根據(jù)反彈的規(guī)則，，在 LB模型中Dt＝1時(shí)，固壁邊界點(diǎn)所受作用力大小為．

1.3 大渦模擬在LBM中的實(shí)現(xiàn)

大渦模擬(large eddy simulation，LES)方法是解決高雷諾數(shù)流動(dòng)模擬的常用方法之一，其主要思想是直接模擬大尺度的紊動(dòng)，而小尺度紊動(dòng)的影響通過(guò)亞格子模型反映，既可以得到比雷諾平均方法更多的大尺度渦結(jié)構(gòu)的信息，又可以比直接數(shù)值模擬節(jié)省大量的計(jì)算量[26]．在 LES模型中通過(guò)過(guò)濾尺度區(qū)分直接求解的大渦尺度和待?；男u尺度，在LB模型中過(guò)濾尺度常選為格子長(zhǎng)度．LB模型中實(shí)現(xiàn)大渦模擬需要引入紊動(dòng)黏度，以反映小尺度渦對(duì)大尺度渦的影響，并通過(guò)將松弛時(shí)間變?yōu)榭紤]了紊動(dòng)黏度影響的總松弛時(shí)間以增大模型中總的黏度.總的黏度由和相加而得

紊動(dòng)黏度通常由速度特征尺度與長(zhǎng)度特征尺度的乘積確定．假設(shè)在亞格子尺度動(dòng)能方程中能量的產(chǎn)生和耗散作用平衡，速度尺度為，長(zhǎng)度尺度為為L(zhǎng)ES模型算子，紊動(dòng)黏度為

式中CX為 LES模型的常數(shù)．加入了紊動(dòng)黏度之后，總松弛時(shí)間為

將式(10)和(11)代入式(12)中可得

不同的LES模型所選取的模型算子不同，其中最易于實(shí)現(xiàn)的是 Smagorinsky模型．其模型算子由速度張量的對(duì)稱部分．

式中Cs為 Smagorinsky常數(shù)；為可解尺度流動(dòng)的應(yīng)變率張量；為網(wǎng)格濾波算符．

然而 Smagorinsky模型具有一些天然的缺陷，如模型參數(shù)Cs需要先驗(yàn)確定，以及在層流區(qū)紊動(dòng)黏度仍不為0，不符合物理實(shí)際等．

WALE模型是一種新的壁面自適應(yīng)局部渦黏模型，可以更好地反映近壁紊動(dòng)特性[16]．該模型算子在引入速度梯度張量的對(duì)稱部分，即變形率張量的基礎(chǔ)上，又引入了速度梯度張量的不對(duì)稱部分，即轉(zhuǎn)動(dòng)張量，其模型算子為

式中轉(zhuǎn)動(dòng)張量為

模型常數(shù).WALE模型的準(zhǔn)確性和精度已經(jīng)過(guò)大量研究驗(yàn)證[16-17,27]，在基于 TRT模型[22]的計(jì)算程序中分別加入 Smagorinsky模型和 WALE大渦模擬模型，對(duì)沙紋床面振蕩流邊界層的紊動(dòng)流場(chǎng)進(jìn)行模擬，對(duì)結(jié)果進(jìn)行比較．

2 模型設(shè)置

選取Freds?e等[1]固定沙紋波浪水槽實(shí)驗(yàn)數(shù)據(jù)進(jìn)行模擬．模型參數(shù)按照實(shí)驗(yàn)參數(shù)設(shè)置，以水質(zhì)點(diǎn)軌道運(yùn)動(dòng)幅值a為特征長(zhǎng)度，振蕩流水質(zhì)點(diǎn)運(yùn)動(dòng)速度幅值U0為特征速度，無(wú)量綱化的沙紋長(zhǎng)度/a＝2.2，沙紋坡度為沙紋波高，振幅雷諾數(shù)Rea＝23,163．計(jì)算區(qū)域順流向長(zhǎng)度Lx=4.4a，即區(qū)域內(nèi)包含兩個(gè)完整沙紋，可以保證計(jì)算區(qū)域?qū)挾茸阋园ㄝ^大尺度渦旋結(jié)構(gòu)．計(jì)算區(qū)域頂邊界設(shè)為自由滑移邊界，與順流向垂直的兩個(gè)邊界和兩個(gè)側(cè)邊界均設(shè)置為周期對(duì)稱邊界，底邊界即沙紋床面表面設(shè)置為無(wú)滑移的固壁邊界．沙紋表面由圓弧組成，沙紋頂峰尖銳，谷底平坦．計(jì)算區(qū)域示意圖見(jiàn)圖 2．模型各方向的格子數(shù)Nx、Ny和Nz分別為 432、412和 432．振蕩流流動(dòng)采用周期性的壓力梯度驅(qū)動(dòng)．算例共模擬 25個(gè)振蕩周期，其中前 10個(gè)周期用來(lái)獲得穩(wěn)定的流場(chǎng)，后面15個(gè)周期用于結(jié)果分析．算例在天河1號(hào)超算中心運(yùn)行，使用72個(gè)CPU，每個(gè)算例運(yùn)行約77h．

/＝0.159，h

圖2 計(jì)算區(qū)域示意Fig.2 Computational domain

3 結(jié)果和分析

3.1 平均流速分布

對(duì)瞬時(shí)水平流速的垂向分布進(jìn)行系綜平均和空間平均，得到圖3所示4個(gè)特征位置處的平均流速分布．圖 4顯示了不同沙紋位置處平均流速的垂向分布與實(shí)驗(yàn)結(jié)果的比較，其中模擬值包括WALE模型和Smagorinsky模型的結(jié)果．圖中縱坐標(biāo)的起點(diǎn)為計(jì)算底邊界，即沙紋谷底位置．實(shí)線為WALE模型結(jié)果，虛線為 Smagorinsky模型結(jié)果，圓圈為實(shí)驗(yàn)結(jié)果．受沙紋形狀的影響，在各沙紋位置處的流速分布呈現(xiàn)出明顯的區(qū)別．兩種大渦模擬模型的結(jié)果存在較明顯的差別，Smagorinsky模型在各沙紋表面位置的多個(gè)相位結(jié)果偏大，而 WALE模型的平均流速分布更接近實(shí)驗(yàn)結(jié)果．

圖3 沙紋位置示意Fig.3 Locations of sand ripple

3.2 紊動(dòng)強(qiáng)度

圖5為采用WALE模型和Smagorinsky模型得到的沙紋各位置順流向平均紊動(dòng)強(qiáng)度的垂向分布與實(shí)驗(yàn)結(jié)果的比較．圖中縱坐標(biāo)的起點(diǎn)為計(jì)算底邊界，即沙紋谷底位置．在沙紋 1/4位置，如圖5(a)中°時(shí)刻，Smagorinsky模型結(jié)果中紊動(dòng)強(qiáng)度大值出現(xiàn)的范圍明顯大于實(shí)驗(yàn)結(jié)果，而 WALE模型的結(jié)果與實(shí)驗(yàn)結(jié)果更加接近．在沙紋谷底位置，如圖5(d)中時(shí)刻，紊動(dòng)強(qiáng)度分布具有雙峰的特征[1,9]，其中下方的峰值產(chǎn)生于邊界層，上方的峰值由沙紋頂峰處的剪切層分離造成．WALE模型結(jié)果盡管幅值大于實(shí)驗(yàn)結(jié)果，但很好地反映了該位置處的這一特征，而 Smagorinsky模型結(jié)果不僅明顯偏小，且沒(méi)有能夠反映雙峰特征．簡(jiǎn)單的 Smagorinsky模型在應(yīng)用于沙紋床面的邊界層流動(dòng)模擬時(shí)，沙紋背流側(cè)坡面和谷底位置的平均流速分布和紊動(dòng)強(qiáng)度結(jié)果均與實(shí)驗(yàn)結(jié)果存在較明顯的差異．相比Smagorinsky模型，WALE模型不僅考慮了應(yīng)變率張量，而且包含了轉(zhuǎn)動(dòng)張量的影響，更適用于復(fù)雜邊界流動(dòng)的模擬，較好地反映了邊界層內(nèi)，特別是沙紋背流側(cè)流動(dòng)分離區(qū)域的紊動(dòng)特征．根據(jù) WALE模型和Smagorinsky模型的比較結(jié)果，在沙紋床面振蕩流邊界層的模擬中，WALE模型可以更好地反映沙紋這種尺度更大的粗糙單元上方邊界層的流動(dòng)特性和紊動(dòng)特性．下文結(jié)果來(lái)自于WALE模型的模擬．

圖4 沙紋各位置處平均流速分布Fig.4 Average velocity distributions at different locations of sand ripple

3.3 壓力分布和形狀摩阻系數(shù)

沙紋的存在使邊界層在沙紋頂峰背后產(chǎn)生流動(dòng)分離，形成渦旋，而渦旋在沙紋背流側(cè)的形成和發(fā)展使沙紋背流面受到負(fù)壓作用．沙紋床面所受的形狀拖曳力Fp就是由沙紋兩側(cè)坡面所受的壓力差產(chǎn)生，可以通過(guò)沙紋表面各處正向壓力水平分量的積分來(lái)計(jì)算．無(wú)量綱的形狀摩阻系數(shù)定義為

圖5 沙紋各位置處平均紊動(dòng)強(qiáng)度分布Fig.5 Average turbulent intensity distributions at different locations of sand ripple

式中A為沙紋床面的水平面積．圖 6為形狀摩阻系數(shù)的歷時(shí)曲線，圖7顯示了振蕩周期中各相位時(shí)刻沙紋附近無(wú)量綱壓力的分布云圖．將圖 7中由左向右方向記為振蕩流的正向流動(dòng)方向，由右向左為反向．在自由流向正向轉(zhuǎn)向時(shí)刻，如圖 7(a)，沙紋背流側(cè)上方存在較大的負(fù)壓區(qū)，該負(fù)壓區(qū)是由脫離沙紋表面的渦旋運(yùn)動(dòng)至該位置造成的．由于該渦旋已經(jīng)不與床面接觸，因此基本沒(méi)有對(duì)床面壓力產(chǎn)生影響.在自由流向正向加速階段，如圖 7(b)，沙紋所受正向壓力的范圍從沙紋頂峰到前一個(gè)沙紋背流側(cè)坡面底部，而從頂峰到背流側(cè)坡面中部，流動(dòng)分離產(chǎn)生的渦旋對(duì)該段沙紋表面造成負(fù)壓作用，沙紋所受壓力的水平合力，即形狀拖曳力符號(hào)為正．在自由流流速達(dá)到幅值的時(shí)刻，如圖 7(c)，沙紋迎流側(cè)坡面所受正壓的范圍和大小均有所減少．在自由流速減速時(shí)刻，如圖 7(d)，受渦旋影響的負(fù)壓區(qū)已經(jīng)脫離沙紋表面，沒(méi)有對(duì)床面產(chǎn)生作用，同時(shí)由于底部流速已經(jīng)提前于自由流流速接近為 0，在沙紋迎流側(cè)坡面沒(méi)有造成正壓作用，因此tw=135°時(shí)刻沙紋兩側(cè)壓力差接近為0．由于流動(dòng)分離、渦旋發(fā)展，導(dǎo)致沙紋兩側(cè)存在隨相位變化的壓力差，這是產(chǎn)生床面形狀拖曳力的直接原因，也是影響床面總拖曳力的重要因素之一．

圖6 形狀摩阻系數(shù)歷時(shí)曲線Fig.6 Time serials of form friction factor

圖7 各時(shí)刻壓力場(chǎng)云圖Fig.7 Pressure contour at different time series

3.4 沙紋床面受力

沙紋床面受到水平方向的總摩阻力FD包括兩部分作用，分別是膚面摩阻力fF和形狀摩阻力Fp．其中膚面摩阻力fF為沙紋表面膚面摩阻的水平方向分量的積分，常用無(wú)量綱的膚面摩阻系數(shù)衡量fF的大小，表達(dá)式為

圖8 膚面摩阻系數(shù)歷時(shí)曲線Fig.8 Time serials of skin friction factor

圖 8為膚面摩阻系數(shù)的歷時(shí)曲線．表1中列出了本節(jié)驗(yàn)證算例和 Dimas等[28]算例的參數(shù)和作用力系數(shù)的結(jié)果．需要說(shuō)明的是，DK算例沙紋坡度為0.14°，略小于本文的 0.159．從結(jié)果可見(jiàn)，本文與 DK的膚面摩阻系數(shù)ff分別為 0.04和 0.03，較為接近．本文形狀摩阻系數(shù)的最大值fp,m為 0.39，略大于DK結(jié)果 0.34，可能的原因是本文算例的沙紋坡度為0.159，大于 DK 算例的 0.140，文獻(xiàn)曾指出形狀摩阻會(huì)隨沙紋坡度增大而增大[29]，因此結(jié)果較合理．本文與 DK算例形狀摩阻系數(shù)與膚面摩阻系數(shù)之比分別為9.8和11.3，證明在該條件下形狀摩阻的作用遠(yuǎn)大于膚面摩阻，與文獻(xiàn)中的結(jié)論相符[30]．由上述比較可知，通過(guò)本文的三維LB模擬可以得到比較合理的沙紋床面形狀摩阻和膚面摩阻結(jié)果．

表1 摩阻系數(shù)與前人數(shù)模結(jié)果的比較Tab.1 Fiction factor result compared with numerical simulation result from literature

3.5 三維渦動(dòng)力特性

圖 9展示了振蕩周期內(nèi)時(shí)刻沙紋上的渦結(jié)構(gòu)．圖中的渦旋結(jié)構(gòu)采用等值面識(shí)別[31]．在時(shí)刻，即圖 9(a)、(b)所示，自由流流速處于正向加速階段，沙紋頂峰背后已經(jīng)發(fā)生流動(dòng)分離，在頂峰背后形成沿z軸方向接近二維的滾動(dòng)渦旋，并向下游運(yùn)動(dòng)．在時(shí)刻，即圖 9(c)、(d)所示，自由流流速達(dá)到幅值，床面上的渦旋結(jié)構(gòu)經(jīng)過(guò)相互作用和與床面的作用，橫向滾動(dòng)渦已經(jīng)分裂成大量小尺度的渦旋結(jié)構(gòu)，在沙紋谷底上方，存在尺度相對(duì)較大的縱向渦結(jié)構(gòu)，而此時(shí)渦旋還沒(méi)有運(yùn)動(dòng)至下一個(gè)沙紋的迎流側(cè)坡面．橫向滾動(dòng)渦出現(xiàn)在沙紋頂峰之后，沿著背浪面發(fā)展，一般存在于周期中加速相位的前期．橫向滾動(dòng)渦旋的出現(xiàn)，使靠近沙紋頂峰位置的床面切應(yīng)力增大．圖 10所示為在沙紋頂峰下游方向位置膚面摩阻系數(shù)的歷時(shí)曲線．在時(shí)刻出現(xiàn)切應(yīng)力幅值，與橫向滾動(dòng)渦出現(xiàn)的空間和時(shí)間對(duì)應(yīng)．縱向渦主要出現(xiàn)在沙紋背流側(cè)坡面和谷底上方，一般存在于加速階段的后期一直到流動(dòng)轉(zhuǎn)向．圖11所示為距離沙紋頂峰位置處膚面摩阻系數(shù)的歷時(shí)曲線．可以看到在區(qū)間膚面摩阻系數(shù)從0.03變到－0.05，負(fù)向的膚面摩阻系數(shù)幅值一直持續(xù)到左右時(shí)刻，該段變化與縱向渦出現(xiàn)的空間和時(shí)間對(duì)應(yīng)．

圖9 沙紋床面上的典型渦結(jié)構(gòu)Fig.9 Typical vortex structure over rippled bed

圖10 /12位置膚面摩阻系數(shù)歷時(shí)曲線Fig.10 Time series of skin friction factor at/12

圖11 /3位置膚面摩阻系數(shù)歷時(shí)曲線Fig.11 Time series of skin friction factor at/3

4 結(jié) 語(yǔ)

將 WALE模型引入三維 LB模型，對(duì)沙紋床面上的振蕩紊流流邊界層進(jìn)行了模擬，得到了與實(shí)驗(yàn)符合較好的水動(dòng)力特性結(jié)果．比較Smagorinsky模型和WALE模型的模擬結(jié)果可知，在沙紋床面的模擬中，WALE模型可以更好地反映振蕩流作用下的紊動(dòng)強(qiáng)度等粗糙紊流特性．WALE模型與LB模型相結(jié)合計(jì)算得到的床面形狀摩阻系數(shù)和膚面摩阻系數(shù)結(jié)果也與前人數(shù)模結(jié)果符合較好．計(jì)算結(jié)果還表明，在沙紋床面的情況下，床面總摩阻力中形狀摩阻相對(duì)膚面摩阻占主要部分，與以往的結(jié)論一致．三維模擬還得到了沙紋上不同時(shí)刻出現(xiàn)的橫向滾動(dòng)渦和縱向渦結(jié)構(gòu)及與膚面摩阻系數(shù)變化的關(guān)系．以上結(jié)果表明，所建立的結(jié)合 WALE大渦模擬模型的三維LB模型可以獲得合理的沙紋床面振蕩流邊界層水動(dòng)力和床面受力結(jié)果，可以利用該模型對(duì)沙紋床面振蕩流粗糙紊流邊界層運(yùn)動(dòng)特性進(jìn)行系統(tǒng)研究，并進(jìn)一步獲得有關(guān)沙紋表觀粗糙度和摩阻系數(shù)等的變化規(guī),．

[1] Freds?e J，Andersen K H，Sumer B M.Wave plus current over a ripple-covered bed[J].Coastal Engineering，1999，38(4)：177-221.

[2] Hare J，Hay A E，Zedel L，et al.Observations of the space-time structure of flow，turbulence，and stress over orbital-scale ripples[J].Journal of Geophysical Research，2014，119(3)：1876-1898.

[3] van der Werf J J，Magar V，Malarkey J，et al.2DV modelling of sediment transport processes over full-scale ripples in regular asymmetric oscillatory flow[J].Continental Shelf Research，2008，28(8)：1040-1056.

[4] 蔣昌波，白玉川，趙子丹，等.波浪作用下沙紋床面底層流動(dòng)特性研究[J].水科學(xué)進(jìn)展，2003，14(3)：333-340.

Jiang Changbo，Bai Yuchuan，Zhao Zidan，et al.Study on wave bottom boundary layers over a rippled bed[J].Advances in Water Science，2003，14(3)：333-340(in Chinese).

[5] 蔣昌波，白玉川，趙子丹，等.波浪作用下渦動(dòng)沙紋床面的懸沙運(yùn)動(dòng)數(shù)值研究[J].水利學(xué)報(bào)，2003(3)：93-97.

Jiang Changbo，Bai Yuchuan，Zhao Zidan，et al.Numerical study on suspended sediment movement in wave boundary layers over rippled bed[J].Journal of Hydraulic Engineering，2003(3)：93-97(in Chinese).

[6] 蔣昌波，白玉川，趙子丹，等.沙紋床面上波流共同作用的數(shù)值模擬[J].水利學(xué)報(bào)，2005，36(1)：62-68.

Jiang Changbo，Bai Yuchuan，Zhao Zidan，et al.Numerical simulation for joint action of wave and current on rippled sea bed[J].Journal of Hydraulic Engineering，2005，36(1)：62-68(in Chinese).

[7] 程永舟，王永學(xué)，蔣昌波.非線性波作用下非對(duì)稱沙紋床面流場(chǎng)特性數(shù)值分析[J].大連理工大學(xué)學(xué)報(bào)，2008，48(3)：423-429.

Cheng Yongzhou，Wang Yongxue，Jiang Changbo.Numerical analysis of characteristics of fluid field over asymmetric rippled bed under nonlinear wave[J].Journal of Dalian University of Technology，2008，48(3)：423-429(in Chinese).

[8] Barr B C，Slinn D N，Pierro T，et al.Numerical simulation of turbulent，oscillatory flow over sand ripples[J].Journal of Geophysical Research，2004，109：C09009-1-19.

[9] Grigoriadis D G E，Balaras E，Dimas A A.Coherent structures in oscillating turbulent boundary layers over a fixed rippled bed[J].Flow，Turbulence and Combustion，2013，91(3)：565-585.

[10] Zhang Q H，Sun Y B，Zhang J F.Simulation of oscillatory laminar boundary layer flow via the lattice Boltzmann method[C]// Proceedings of 3rd Internatioal Conference on Asian and Pacific Coasts.Jeju，Korea，2005：271-274.

[11] 孫亞斌，張慶河，張金鳳.振蕩層流邊界層運(yùn)動(dòng)的格子Boltzmann模擬[J].水動(dòng)力學(xué)研究與進(jìn)展，2006，21(3)：347-353.

Sun Yabin，Zhang Qinghe，Zhang Jinfeng.Simulation of oscillatory laminar boundary layer flow based on lattice Boltzmann method[J].Journal of Hydrodynamics，2006，21(3)：347-353(in Chinese).

[12] Cappietti L，Chopard B.A lattice Boltzmann study of the 2D boundary layer created by an oscillating plate[J].International Journal of Modern Physics C，2006，17(1)：39-52.

[13] 丁 磊，張慶河.振蕩流邊界層中顆粒受力的三維格子玻耳茲曼模擬[J].水動(dòng)力學(xué)研究與進(jìn)展，2010，25(3)：391-397.

Ding Lei，Zhang Qinghe.3D lattice Boltzmann simulation of forces on a fixed spherical[J].Journal of Hydrodynamics，2010，25(3)：391-397(in Chinese).

[14] Ding L，Zhang Q H.Lattice Boltzmann simulation to characterize roughness effects of oscillatory boundary layer flow over a rough bed[C]//Proceedings of 32nd Conference on Coastal Engineering.Shanghai，China，2010：1397-1407.

[15] 周志博，張慶河.粗糙床面振蕩紊流邊界層運(yùn)動(dòng)的三維格子玻爾茲曼模擬[J].水動(dòng)力學(xué)研究與進(jìn)展，2016，31(4)：463-471.

Zhou Zhibo，Zhang Qinghe.3-D lattice Boltzmann simulation of turbulent oscillatory boundary layer flow over rough beds[J].Chinese Journal of Hydrodynamics，2016，31(4)：463-471(in Chinese).

[16] Nicoud F，Ducros F.Subgrid-scale stress modelling based on the square of the velocity gradient tensor[J].Flow，Turbulence and Combustion，1999，62(3)：183-200.

[17] Weickert M，Teike G，Schmidt O，et al.Investigation of the LES WALE turbulence model within the lattice Boltzmann framework[J].Computers & Mathematics with Applications，2010，59(7)：2200-2214.

[18] Qian Y H，d'Humières D，Lallemand P.Lattice BGK models for navier-stokes equation[J].Europhysics Letters，1992，17(6)：479-484.

[19] Bhatnagar P L，Gross E P，Krook M.A model for collision processes in gases[J].Physical Review，1954，94：511-524.

[20] d'Humieres D.Generalized lattice-Boltzmann equations[C]//Rarefied Gas Dynamics：Theory and Simulations.Washington，USA，1992：450-458.

[21] Lallemand P，Luo L S.Theory of the lattice Boltzmann method：Dispersion，dissipation，isotropy，Galilean invariance，and stability[J].Physical Review E，2000，61(6)：6546-6562.

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

[23] He X，Zou Q，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.

[24] Nguyen N Q，Ladd A J C.Sedimentation of hard-sphere suspensions at low Reynolds number[J].Journal of Fluid Mechanics，2005，525：73-104.

[25] Ladd A J C.Numerical simulations of particulate suspensions via a discretized Boltzmann equation(Part 1)：Theoretical foundation[J].Journal of Fluid Mechanics，1994，271：285-309.

[26] 崔桂香，許春曉，張兆順.湍流大渦數(shù)值模擬進(jìn)展[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2004，22(2)：121-129.

Cui Guixiang，Xu Chunxiao，Zhang Zhaoshun.Progress in large eddy simulation of turbulent flows[J].Acta Aerodynamica Sinica，2004，22(2)：121-129(in Chinese).

[27] Temmerman L，Leschziner M A，Mellen C P，et al.Investigation of wall-function approximations and subgrid-scale models in large eddy simulation of separated flow in a channel with streamwise periodic constrictions[J].International Journal of Heat and Fluid Flow，2003，24(2)：157-180.

[28] Dimas A A，Kolokythas G A.Flow dynamics and bed resistance of wave propagation over bed ripples[J].Journal of Waterway，Port，Coastal，and Ocean Engineering，2011，137(2)：64-74.

[29] Kim H.Effective form roughness of ripples for waves[J].Journal of Coastal Research，2004，20(3)：731-738.

[30] Wikramanayake P N，Madsen O S.Calculation of Movable Bed Friction Factors[M].McNair E C，Russell K T.Dredging Research Program.Vicksburg：Army Corps of Engineers：1994：1-104.

[31] Jeong J，Hussain F.On the identification of a vortex[J].Journal of Fluid Mechanics，1995，285：69-94.