王 碩, 凡鳳仙, 張盟盟, 胡曉紅

(1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院,上海 200093;2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室,上海 200093)

基于OpenFOAM的正弦形沙丘表面流場特性研究

王 碩1,2, 凡鳳仙1,2, 張盟盟1,2, 胡曉紅1,2

(1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院,上海 200093;2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室,上海 200093)

基于開源計(jì)算流體力學(xué)軟件OpenFOAM,對(duì)沙丘表面流場進(jìn)行數(shù)值模擬,將所得到的沙丘表面風(fēng)速分布與野外觀測結(jié)果進(jìn)行對(duì)比,驗(yàn)證數(shù)值模擬的準(zhǔn)確性.在此基礎(chǔ)上,獲得了距沙丘表面不同高度的速度分布曲線和沙丘表面上的速度矢量場,探討了數(shù)值模擬結(jié)果對(duì)計(jì)算網(wǎng)格的敏感性,考察了沙丘寬高比對(duì)流場特性的影響.結(jié)果表明：OpenFOAM的標(biāo)準(zhǔn)k-ε模型能夠預(yù)測低矮正弦形沙丘表面流場的速度轉(zhuǎn)變特性;在貼近沙丘表面的內(nèi)層,沙丘表面特定高度處的速度分布曲線相對(duì)于沙丘輪廓存在相位提前,隨著高度的增加,相位提前量減少,速度變化趨于平緩;網(wǎng)格數(shù)目對(duì)外層速度分布影響較小,但對(duì)內(nèi)層的速度分布則有重要影響;足夠多的網(wǎng)格數(shù)目是獲得精確的沙丘表面流場信息的必要條件,在沙丘流場數(shù)值模擬中應(yīng)當(dāng)引起重視;保持沙丘寬度不變,隨著沙丘寬高比的減小,沙丘表面速度分布曲線的相位提前量變大,背風(fēng)坡出現(xiàn)回流區(qū),寬高比越小,回流區(qū)的長度與強(qiáng)度越大.

計(jì)算流體力學(xué); 正弦形沙丘; 湍流; 速度轉(zhuǎn)變; 網(wǎng)格敏感性

沙丘是指由于沙粒在流場中相互作用所形成的堆積體,常見于海岸、河谷和干旱沙漠地帶.沙丘表面流場決定著局地沙粒輸運(yùn)的強(qiáng)度,從而影響沙丘的形態(tài),在沙丘演變動(dòng)力學(xué)研究中占據(jù)重要地位,因而備受研究者關(guān)注[1-3].目前對(duì)沙丘表面流場進(jìn)行研究的方法可概括為3大類：野外觀測[4-6]、風(fēng)洞實(shí)驗(yàn)[7-8]和數(shù)值模擬[9-12].野外觀測是最為直接的方法,但是在進(jìn)行大空間尺度和時(shí)間尺度的沙丘研究時(shí),實(shí)驗(yàn)測試和數(shù)據(jù)處理工作量將會(huì)非常復(fù)雜龐大,且難以給出流場結(jié)構(gòu)的詳細(xì)信息.風(fēng)洞實(shí)驗(yàn)可以較好地反映沙丘環(huán)境,但是在滿足相似性條件時(shí),難以同時(shí)保證幾何相似和動(dòng)力學(xué)相似.隨著計(jì)算流體力學(xué)的發(fā)展和日趨成熟,數(shù)值模擬展示出其方便、靈活的優(yōu)勢(shì),不僅能夠?qū)崿F(xiàn)長時(shí)間、大范圍的沙丘流場研究,還可以靈活調(diào)整計(jì)算參數(shù),降低研究的時(shí)間代價(jià)和成本,尤其是在沙丘流場細(xì)節(jié)信息的獲取上,數(shù)值模擬成為一種有效途徑.

1 國內(nèi)外研究現(xiàn)狀

近年來,一些學(xué)者針對(duì)沙丘表面流場特性開展了一系列研究.2004年,Parsons等[13]借助商用軟件PHOENICS,利用RNGk-ε湍流模型對(duì)剖面為三角形的風(fēng)成沙丘表面流場特性進(jìn)行了二維數(shù)值模擬,準(zhǔn)確地預(yù)測到沙丘表面流型,獲得了氣流在坡腳的停滯區(qū)、迎風(fēng)坡的加速區(qū)和背風(fēng)坡的流動(dòng)分離和回流區(qū).2005年,Herrmann等[14]采用商用軟件FLUENT的k-ε模型對(duì)三維新月形沙丘和剖面近似為三角形的二維橫向沙丘表面流場開展數(shù)值模擬,展示了沙丘背風(fēng)坡流動(dòng)分離的現(xiàn)象.2010年,江麗娟等[9]基于FLUENT的k-ε模型對(duì)新月形沙丘表面流場進(jìn)行數(shù)值模擬,所得結(jié)果能夠反映剪應(yīng)力的定性特征,但在定量上與實(shí)驗(yàn)有一定差異.2012年,Balogh等[15]對(duì)湍流模型和壁面函數(shù)進(jìn)行修正,分別基于FLUENT和開源軟件OpenFOAM研究了復(fù)雜地形大氣邊界層流場特性,結(jié)果顯示OpenFOAM具有更高的準(zhǔn)確性和靈活性.2013年,周曉斯等[10]基于大渦模擬(LES)方法,自行編制計(jì)算程序,對(duì)小尺度新月形沙丘背風(fēng)坡流場特性進(jìn)行二維數(shù)值模擬,所得風(fēng)速與風(fēng)洞實(shí)驗(yàn)結(jié)果吻合良好.同年,Araújo等[11]采用FLUENT的k-ε模型,對(duì)中剖面近似為三角形的橫向沙丘表面流場開展數(shù)值模擬,發(fā)現(xiàn)了回流區(qū)長度隨剪切速度的變化關(guān)系.2014年,Peralta等[16]基于OpenFOAM的標(biāo)準(zhǔn)k-ε、RNGk-ε和k-ωSST模型,采用結(jié)構(gòu)化和非結(jié)構(gòu)化兩種網(wǎng)格類型,模擬了兩種復(fù)雜地形上大氣邊界層風(fēng)場,并將時(shí)均風(fēng)速和湍動(dòng)能的模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比,結(jié)果顯示:第一種地形模擬結(jié)果和實(shí)驗(yàn)結(jié)果符合較好,在采用RNGk-ε湍流模型和結(jié)構(gòu)化網(wǎng)格時(shí)吻合最好;而對(duì)于第二種地形,模擬結(jié)果與實(shí)驗(yàn)結(jié)果偏差較大.2015年,Michelsen等[12]基于OpenFOAM的RNGk-ε湍流模型,對(duì)新月形沙丘表面流場開展二維和三維模擬,結(jié)果顯示,二維模擬低估了沙丘表面風(fēng)速.

上述工作主要通過與野外觀測和風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)的對(duì)比來驗(yàn)證模型,進(jìn)而對(duì)沙丘流場動(dòng)力學(xué)開展模擬.然而,不同的研究者所針對(duì)的地形和空間尺度差別很大,采用的CFD軟件及湍流模型也有所不同,所得的結(jié)論缺乏普適性.因此,有必要結(jié)合實(shí)驗(yàn)結(jié)果,對(duì)沙丘表面流場特性開展更細(xì)致的研究,以尋求適應(yīng)性強(qiáng)、快速高效的數(shù)值模擬方法,進(jìn)而獲得準(zhǔn)確的沙丘表面流場信息.OpenFOAM因其免費(fèi)、開源、穩(wěn)定強(qiáng)大的底層類庫,以及能夠并行計(jì)算的優(yōu)點(diǎn)受到CFD用戶的關(guān)注,并在大氣邊界層模擬中發(fā)揮著日益重要的作用[12,15-16].最近,Claudin等[4]首次報(bào)道了低矮正弦形沙丘表面流場的野外觀測結(jié)果,獲得了貼近沙丘表面的風(fēng)速分布,為沙丘表面流場的精細(xì)數(shù)值模擬提供了新的實(shí)驗(yàn)依據(jù).本文借助OpenFOAM的標(biāo)準(zhǔn)k-ε模型,采用與文獻(xiàn)[4]相同的沙丘輪廓,對(duì)沙丘表面流場開展數(shù)值模擬,將數(shù)值模擬結(jié)果與文獻(xiàn)[4]數(shù)據(jù)相結(jié)合,分析流場特性.在此基礎(chǔ)上,探討了網(wǎng)格數(shù)目對(duì)計(jì)算結(jié)果的影響,并對(duì)不同寬高比沙丘表面流場進(jìn)行預(yù)測.本文研究可為進(jìn)一步利用數(shù)值模擬方法研究沙丘表面沙粒輸運(yùn)提供基礎(chǔ)和參考.

2 計(jì)算模型

沙丘表面流場為處于大氣邊界層內(nèi)的近地面流場,其特點(diǎn)是小流速、大雷諾數(shù)的湍流.這種湍流速度場往往具有擾動(dòng)尺度小且變化頻率高的特點(diǎn),在建立描述流場的納維-斯托克斯方程時(shí)通常采用將速度分解為時(shí)均速度和脈動(dòng)速度,以降低求解的時(shí)間成本.將風(fēng)速方向設(shè)定為x向,豎直向上方向設(shè)定為z向,與x和z垂直的方向設(shè)定為y向,則在笛卡爾坐標(biāo)系下雷諾平均納維-斯托克斯(RANS)方程可寫為

式中:t為時(shí)間;u為流速;p為壓力;ρ為空氣密度;μ為空氣動(dòng)力粘度;δij為克羅內(nèi)克函數(shù),當(dāng)i=j時(shí),δij=1,當(dāng)i≠j時(shí),δij=0.

為使方程封閉,引入標(biāo)準(zhǔn)k-ε湍流模型計(jì)算湍動(dòng)能k和耗散率ε,即

其中,

(5)

式中:σk,σε,C1ε,C2ε,Cμ為模型常數(shù),在大氣邊界層模擬中通常采用σk=1.0,σε=1.30,C1ε=1.21,C2ε=1.92,Cμ=0.03[17-19];Eij為平均應(yīng)變率,其表達(dá)式為

(6)

3 計(jì)算條件及方法

3.1 計(jì)算區(qū)域

為了驗(yàn)證計(jì)算模型和計(jì)算方法的準(zhǔn)確性,對(duì)文獻(xiàn)[4]中沙丘表面流場進(jìn)行數(shù)值模擬,并與其野外觀測數(shù)據(jù)進(jìn)行對(duì)比.文獻(xiàn)[4]所研究的沙丘位于平坦的地面上,沙丘輪廓可以近似為正弦曲線

zd=zref+Acos(Kx)

(7)

式中:zd為沙丘表面的z坐標(biāo);zref為參考高度;A為幅值;K為波數(shù),K=2π/λ,λ為波長;x=0對(duì)應(yīng)于沙丘坡頂位置,zd=0對(duì)應(yīng)于沙丘坡腳位置.文獻(xiàn)[4]給出,zref=A≈1.1 m,K≈0.14 m-1,λ≈46 m.此外,為通過數(shù)值模擬考察沙丘寬高比對(duì)流場特性的影響,計(jì)算時(shí)保持沙丘寬度(即波長λ)不變,通過調(diào)整沙丘高度(2A)以改變沙丘的寬高比.

開展數(shù)值模擬首先應(yīng)當(dāng)確定計(jì)算區(qū)域,并劃分計(jì)算網(wǎng)格.由于沙丘表面速度變化集中在來流方向(x向)和沙丘高度方向(z向),沙丘長度方向(y向)速度變化很小;為減少網(wǎng)格數(shù)目,縮短計(jì)算時(shí)間,開展二維模擬.為避免邊界效應(yīng),進(jìn)口、出口、頂面與沙丘的距離應(yīng)足夠大.為此,選擇計(jì)算區(qū)域?qū)挾圈=200 m,計(jì)算區(qū)域高度根據(jù)沙丘坡頂高度進(jìn)行調(diào)整,保持兩高度之比不發(fā)生變化,其中沙丘幅值A(chǔ)=1.1 m時(shí),Δz=20 m,沙丘位于計(jì)算區(qū)域中心,如圖1所示.

圖1 計(jì)算區(qū)域示意圖Fig.1 Schematic of the simulation domain

3.2 網(wǎng)格劃分

計(jì)算網(wǎng)格由OpenFOAM中的blockMesh和snappyHexMesh網(wǎng)格生成程序產(chǎn)生.在數(shù)值模擬中,網(wǎng)格越小,計(jì)算精度越高,但計(jì)算耗時(shí)也越長.為解決計(jì)算精度和計(jì)算速度之間的矛盾,本文分區(qū)域設(shè)置網(wǎng)格大小等級(jí),對(duì)于速度變化較大的沙丘附近區(qū)域采用較密的網(wǎng)格.以幅值A(chǔ)=1.1 m的沙丘表面流場計(jì)算區(qū)域?yàn)槔?圖2給出了不同等級(jí)網(wǎng)格區(qū)域分布,在研究不同高度沙丘時(shí),各區(qū)域高度按比例調(diào)整.總的來說,區(qū)域I網(wǎng)格最大,區(qū)域Ⅱ網(wǎng)格尺寸(長和高)為區(qū)域I的1/2,區(qū)域Ⅲ網(wǎng)格尺寸為區(qū)域Ⅱ的1/2.在計(jì)算區(qū)域底面(即包含沙丘的地面)添加5層邊界層,在遠(yuǎn)離下表面的方向上,邊界層網(wǎng)格厚度呈比例增加,比例系數(shù)為1.1.

圖2 不同等級(jí)網(wǎng)格區(qū)域分布Fig.2 Distribution of regions with different mesh levels

3.3 邊界條件

進(jìn)口采用速度邊界條件,速度值可由大氣邊界層風(fēng)速的對(duì)數(shù)分布率給出[4]

(8)

式中:κ為馮·卡門常數(shù),κ=0.41;u*為剪切速度,u*=0.5 m/s;z0為空氣動(dòng)力學(xué)粗糙長度,z0=0.83 mm.

湍動(dòng)能和耗散率的表達(dá)式為[20-21]

(9)

(10)

底面采用壁面邊界條件;頂面設(shè)置為移動(dòng)壁面,頂面風(fēng)速由式(8)確定;出口設(shè)置為壓力出口.此外,對(duì)于本文的二維數(shù)值模擬,需要在y向采用單層網(wǎng)格,并將相應(yīng)的邊界設(shè)置為“empty”.

3.4 求解控制

基于SIMPLE算法[22-24]對(duì)壓力和動(dòng)量方程進(jìn)行耦合求解,該算法利用有限體積方法求解RANS方程,由于在計(jì)算過程中標(biāo)量儲(chǔ)存在網(wǎng)格中心,矢量儲(chǔ)存在網(wǎng)格面中心,因而需要采用插值方法.本文的數(shù)值計(jì)算中,梯度項(xiàng)、拉普拉斯項(xiàng)采用高斯線性插值格式;速度項(xiàng)采用高斯QUICKV插值格式;湍動(dòng)能項(xiàng)以及湍動(dòng)能擴(kuò)散項(xiàng)則采用高斯迎風(fēng)格式.壓力、速度、湍動(dòng)能和耗散率的松弛因子分別為0.3,0.7,0.7和0.7.收斂條件為質(zhì)量和動(dòng)量通量的殘差達(dá)到10-8,壓力的殘差達(dá)到10-7.

4 結(jié)果與討論

4.1 數(shù)值模擬結(jié)果的驗(yàn)證

為驗(yàn)證數(shù)學(xué)模型和數(shù)值計(jì)算方法的合理性,將數(shù)值模擬得到的沙丘表面速度分布與文獻(xiàn)[4]的野外觀測結(jié)果作對(duì)比,如圖3所示.其中,計(jì)算區(qū)域網(wǎng)格數(shù)目為1.34×105個(gè),圖3(a)給出了距離沙丘表面0.11 m高度處(內(nèi)層)的速度分布,圖3(b)給出了距離沙丘表面0.50 m高度處(外層)的速度分布.圖中,u0為最大速度和最小速度的平均值,連續(xù)的曲線為數(shù)值模擬結(jié)果,帶誤差條的離散點(diǎn)為來自文獻(xiàn)[4]的野外觀測結(jié)果.可見,數(shù)值模擬和文獻(xiàn)[4]結(jié)果吻合較好,表明二維模擬能夠預(yù)測正弦形沙丘表面流場特性.但也需要注意到,在一些位置模擬值比文獻(xiàn)[4]略低,這可能是由于模擬時(shí)將三維沙丘視為二維所致,已有對(duì)新月形沙丘表面流場的數(shù)值模擬結(jié)果顯示,二維模擬低估了沙丘表面風(fēng)速[12],后續(xù)研究可采用三維模擬進(jìn)一步優(yōu)化計(jì)算結(jié)果.圖3的數(shù)值模擬和文獻(xiàn)[4]結(jié)果均反映出:內(nèi)層氣流在迎風(fēng)坡加速,在坡頂上游速度達(dá)到最大值,速度分布相對(duì)沙丘輪廓具有相位提前量φ/K≈3.2 m;外層氣流速度相位提前量明顯減小.內(nèi)層和外層速度分布差異是由于流體力學(xué)機(jī)理不同所致[25-27].在內(nèi)層,流體粘性起重要作用,由于粘性阻力的影響,風(fēng)速峰值出現(xiàn)在沙丘坡頂前;隨著高度的增加,慣性作用增強(qiáng),粘性的影響減弱,風(fēng)速曲線的相位提前量減小.沙丘表面流場內(nèi)層速度具有相位提前的現(xiàn)象,這對(duì)沙丘的形成、移動(dòng)以及演變具有重要影響.氣流在迎風(fēng)側(cè)加速,背風(fēng)側(cè)減速,造成表面剪切應(yīng)力在迎風(fēng)側(cè)增加,背風(fēng)側(cè)減小,從而沙丘上游產(chǎn)生吹蝕區(qū),下游產(chǎn)生沉積區(qū),進(jìn)而引發(fā)沙丘的演變.

4.2 沙丘表面流場速度分布

圖4給出了沙丘表面0.11,0.50,1.50,2.00 m高度處的速度分布,數(shù)值模擬采用的網(wǎng)格數(shù)目與圖3相同.可以看出,速度分布隨高度發(fā)生顯著變化,隨著高度的增加,最大速度出現(xiàn)的位置向坡頂遷移,沙丘迎風(fēng)坡和背風(fēng)坡速度趨于對(duì)稱分布,速度變化變得平緩.

圖5給出了沙丘表面流場速度矢量圖.由圖可見,受沙丘表面摩擦的影響,貼近沙丘表面的區(qū)域速度較低,且在迎風(fēng)坡具有向上的分量,在背風(fēng)坡具有向下的分量;在離沙丘表面較遠(yuǎn)的區(qū)域,沙丘表面的影響減弱,流場速度的豎直分量基本消失,速度方向與來流方向平行.已有研究曾對(duì)沙丘背風(fēng)坡流動(dòng)分離的現(xiàn)象作過報(bào)道[11,13-14,28-29],流動(dòng)分離的形成與沙丘大小和外形、來流速度均有關(guān)系.在本節(jié)研究所針對(duì)的低矮沙丘及中等來流速度條件下,無論野外觀測還是數(shù)值模擬均未發(fā)現(xiàn)流動(dòng)分離現(xiàn)象.由圖5還可以看出,數(shù)值模擬可以全面展現(xiàn)流場豐富的細(xì)節(jié)信息,是對(duì)野外觀測的有效補(bǔ)充.

圖4 沙丘表面不同高度處的速度分布Fig.4 Velocity profiles at different heights over the dune

4.3 網(wǎng)格數(shù)目對(duì)流場的影響

為考察數(shù)值模擬結(jié)果對(duì)計(jì)算網(wǎng)格的敏感性,圖6給出了計(jì)算區(qū)域總網(wǎng)格數(shù)目不同時(shí)得到的沙丘表面流場速度分布.其中,圖6(a),(b)分別為沙丘表面0.11 ,0.50 m處的結(jié)果.對(duì)比可知,網(wǎng)格數(shù)目對(duì)內(nèi)層流場速度有重要影響,而對(duì)外層流場速度的影響較小.內(nèi)層速度分布曲線相對(duì)沙丘輪廓曲線具有相位提前量的特征已被野外觀測所證實(shí)[4],并能通過理論分析得到合理解釋[30].然而,在網(wǎng)格數(shù)目較少,例如為4.10×104,7.30×104個(gè)時(shí),數(shù)值模擬結(jié)果無法反映出速度相位提前這一重要特征.外層速度分布曲線雖然對(duì)網(wǎng)格數(shù)目敏感性較弱,但沙丘坡頂附近以及背風(fēng)坡流速隨著網(wǎng)格數(shù)目的增加略有增大.同樣反映出網(wǎng)格數(shù)目較少時(shí),無法體現(xiàn)速度分布的不對(duì)稱性.圖6中4種網(wǎng)格數(shù)目情況下,網(wǎng)格數(shù)目為1.34×105個(gè)時(shí)的數(shù)值模擬結(jié)果與野外觀測結(jié)果吻合良好,且繼續(xù)增大網(wǎng)格數(shù)目對(duì)計(jì)算結(jié)果的影響很小,是數(shù)值模擬中適宜的網(wǎng)格選擇.該網(wǎng)格數(shù)目條件下,區(qū)域Ⅰ,Ⅱ,Ⅲ (見圖2)中網(wǎng)格尺寸分別為0.2 m×0.2 m,0.1 m×0.1 m,0.05 m×0.05 m.本文采用的分區(qū)劃分網(wǎng)格方法,以及網(wǎng)格相對(duì)沙丘寬度和高度的尺寸,可為正弦形表面流場特性數(shù)值模擬研究中的網(wǎng)格劃分提供參考.

圖5 沙丘表面流場速度矢量圖Fig.5 Velocity vectors of the flow field over the dune

圖6 不同網(wǎng)格數(shù)目時(shí)的速度分布Fig.6 Velocity profiles for different mesh numbers

4.4 沙丘寬高比對(duì)流場的影響

為研究正弦形沙丘幾何參數(shù)對(duì)流場特性的影響,對(duì)不同寬高比的沙丘表面流場進(jìn)行數(shù)值預(yù)測.在A=1.1 ,1.65 ,2.2 ,2.75 ,3.3 m,即寬高比λ/(2A)=20.9,13.9,10.5,8.36,6.97情況下獲得了沙丘表面流場信息.圖7給出了距沙丘表面0.11 m與0.50 m處的流場速度分布情況.由圖可知,在距離沙丘表面0.11 m處,隨著沙丘寬高比的減小,一方面,速度分布曲線相對(duì)于沙丘輪廓線的相位提前量增大,迎風(fēng)坡氣流速度變化加劇;另一方面,當(dāng)沙丘寬高比λ/(2A)≤10.5時(shí),背風(fēng)坡氣流速度出現(xiàn)負(fù)值,表明氣流流過坡頂后下沉形成回流區(qū).沙丘迎風(fēng)坡的氣流加速和背風(fēng)坡的回流區(qū)與沙丘侵蝕率有重要關(guān)聯(lián),沙丘寬高比越小,侵蝕越顯著,沙丘移動(dòng)和演變?cè)窖杆?在距離沙丘表面0.5 m處,速度分布曲線與沙丘輪廓曲線的相位差均較小,尤其是在沙丘寬高比λ/(2A)≥10.5時(shí),相位差幾乎為0;同時(shí),沙丘寬高比越小,迎風(fēng)坡氣流加速越快.沙丘寬高比λ/(2A)≥10.5時(shí),在背風(fēng)坡氣流速度全為正;而λ/(2A)≤8.36時(shí),可觀察到速度負(fù)值,且高寬比越小,速度大小越偏離0.這些結(jié)果表明,隨著沙丘寬高比的減小,回流區(qū)長度和強(qiáng)度均增大.不同尺度沙丘表面流場特性研究是后續(xù)利用氣固兩相流數(shù)值模擬方法考察沙粒輸運(yùn),進(jìn)而認(rèn)識(shí)和掌握沙丘的形成機(jī)制、發(fā)展規(guī)律及變化趨勢(shì)的基礎(chǔ).

圖7 不同寬高比沙丘表面的速度分布Fig.7 Velocity profiles over dunes with different aspect ratios

5 結(jié) 論

本文基于開源CFD軟件OpenFOAM對(duì)正弦形沙丘表面流場進(jìn)行了數(shù)值模擬,并將數(shù)值模擬結(jié)果與文獻(xiàn)[4]的野外觀測結(jié)果相對(duì)比,驗(yàn)證數(shù)學(xué)模型和計(jì)算方法的正確性.在此基礎(chǔ)上,對(duì)沙丘表面流場特性進(jìn)行分析,探討了數(shù)值模擬結(jié)果的網(wǎng)格敏感性,考察了沙丘寬高比對(duì)流場特性的影響,得到以下結(jié)論：

a. 基于OpenFOAM的標(biāo)準(zhǔn)k-ε模型能夠得到低矮正弦形沙丘表面流場特性,由其預(yù)測出的內(nèi)層和外層速度分布與文獻(xiàn)[4]中的野外觀測結(jié)果吻合良好;

b. 在貼近沙丘表面的內(nèi)層,速度分布曲線相對(duì)于沙丘輪廓線存在相位提前,隨著高度的增加,相位提前量減少,速度變化趨于平緩,直至風(fēng)速的豎直分量不再存在,沙丘形狀對(duì)流場的影響消失;

c. 計(jì)算網(wǎng)格數(shù)目對(duì)外層速度分布影響較小,但對(duì)內(nèi)層的速度分布則有重要影響,足夠多的網(wǎng)格數(shù)目是獲得精確的沙丘表面流場信息的必要條件,在沙丘流場數(shù)值模擬研究中應(yīng)當(dāng)引起重視;

d. 沙丘寬度不變時(shí),沙丘寬高比減小,沙丘表面速度分布曲線的相位提前量增加,背風(fēng)坡出現(xiàn)回流區(qū),并且寬高比越小,回流區(qū)的長度與強(qiáng)度越大.

[1] PARTELI E J R,KROY K,TSOAR H,et al.Morphodynamic modeling of aeolian dunes:review and future plans[J].The European Physical Journal Special Topics,2014,223(11):2269-2283.

[2] CHARRU F,ANDREOTTI B,CLAUDIN P.Sand ripples and dunes[J].Annual Review of Fluid Mechanics,2013,45(1):469-493.

[3] 呂萍,董治寶.沙丘動(dòng)力學(xué)數(shù)值模擬研究方法綜述(Ⅰ):沙丘周圍流場模擬[J].中國沙漠,2015,35(3):526-533.

[4] CLAUDIN P,WIGGS G F S,ANDREOTTI B.Field evidence for the upwind velocity shift at the crest of low dunes[J].Boundary-Layer Meteorology,2013,148(1):195-206.

[5] PELLETIER J D.Deviations from self-similarity in Barchan form and flux:the case of the Salton sea dunes,California[J].Journal of Geophysical Research:Earth Surface,2013,118(4):2406-2420.

[6] 張正偲,董治寶.橫向沙丘表面氣流湍流特征野外觀測[J].地理科學(xué),2015,35(5):652-657.

[7] CREYSSELS M,DUPONT P,MOCTAR A,et al.Saltating particles in a turbulent boundary layer:experiment and theory[J].Journal of Fluid Mechanics,2009,625(12):47-74.

[8] 賈攀,王元,張洋.縮尺沙丘演化特性的風(fēng)洞實(shí)驗(yàn)研究[J].西安交通大學(xué)學(xué)報(bào),2013,47(5):67-71.

[9] 江麗娟,馬高生.新月形沙丘表面流場的數(shù)值模擬[J].蘭州大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,46(1):48-53.

[10] 周曉斯,王元,李志強(qiáng).小尺度新月形沙丘背風(fēng)側(cè)流場特性的大渦模擬分析[J].西安交通大學(xué)學(xué)報(bào),2013,47(7):114-119.

[12] MICHELSEN B,STROBL S,PARTELI E R J,et al.Two-dimensional airflow modeling under predicts the wind velocity over dunes[J].Scientific Reports,2015,5:16572.

[13] PARSONS D R,WIGGS G F S,WALKER I J,et al.Numerical modelling of airflow over an idealised transverse dune[J].Environmental Modelling & Software,2004,19(2):153-162.

[14] HERRMANN H J,ANDRADE JR J S,SCHATZ V,et al.Calculation of the separation streamlines of barchans and transversedunes[J].Physica A:Statistical Mechanics and Its Applications,2005,357(1):44-49.

[15] BALOGH M,PARENTE A,BENOCCI C.RANS simulation of ABL flow over complex terrains applying an Enhancedk-εmodel and wall function formulation:implementation and comparison for fluent and OpenFOAM[J].Journal of Wind Engineering and Industrial Aerodynamics,2012(104/105/106):360-368.

[16] PERALTA C,NUGUSSE H,KOKILAVANI S P,et al.Validation of thesimpleFoam (RANS) solver for the atmospheric boundary layer in complex terrain[C]//ITM Web of Conferences.London:EDP Sciences,2014,2:01002.

[17] HARGREAVES D M,WRIGHT N G.On the use of thek-εmodel in commercial CFD software to model the neutral atmospheric boundary layer[J].Journal of Wind Engineering and Industrial Aerodynamics,2007,95(5):355-369.

[18] ZHANG X D.CFD simulation of neutral ABL flows[R].Roskilde:Technical University of Denmark,2009.

[19] TAPIA X P.Modelling of wind flow over complex terrain using OpenFoam[D].G?vle:University of G?vle,2009:35-42.

[20] BLOCKEN B,STATHOPOULOS T,CARMELIET J.CFD simulation of the atmospheric boundary layer:wall function problems[J].Atmospheric Environment,2007,41(2):238-252.

[21] YAN B W,LI Q S,HE Y C,et al.RANS simulation of neutral atmospheric boundary layer flows over complex terrain by proper imposition of boundary conditions and modification on thek-εmodel[J].Environmental Fluid Mechanics,2016,16(1):1-23.

[22] VERSTEEG H K,MALALASEKERA W.An introduction to computational fluid dynamics:the finite volumemethod[M].2nd ed.Harlow:Pearson Education Limited,2007.

[23] JASAK H.Error analysis and estimation for the finite volume method with applications to fluid flows[D].London:Imperial College,1996.

[24] 陶文銓.數(shù)值傳熱學(xué)[M].2版.西安:西安交通大學(xué)出版社,2001.

[25] WEAVER C M,WIGGS G F S.Field measurements of mean and turbulent airflow over abarchan sand dune[J].Geomorphology,2011,128(1/2):32-41.

[26] BRITTER R E,HUNT J C R,RICHARDS K J.Air flow over a two-dimensional hill:studies of velocity speed-up,roughness effects andturbulence[J].Quarterly Journal of the Royal Meteorological Society,1981,107(451):91-110.

[27] BADDOCK M C,WIGGS G F S,LIVINGSTONE I.A field study of mean and turbulent flow characteristics upwind,over and downwind of barchan dunes[J].Earth Surface Processes and Landforms,2011,36(11):1435-1448.

[28] WALKER I J,NICKLING W G.Dynamics of secondary airflow and sediment transport over and in the lee of transversedunes[J].Progress in Physical Geography,2002,26(1):47-75.

[29] GROH C,WIERSCHEM A,AKSEL N,et al.Barchan dunes in two dimensions:experimental tests for minimal models[J].Physical Review E,2008,78(2):021304.

[30] COLOMBINI M.Revisiting the linear theory of sand dune formation[J].Journal of Fluid Mechanics,2004,502:1-16.

(編輯:丁紅藝)

Flow Behavior over a Sinusoidal Dune Based on OpenFOAM

WANG Shuo1,2, FAN Fengxian1,2, ZHANG Mengmeng1,2, HU Xiaohong1,2

(1.SchoolofEnergyandPowerEngineering,UniversityofShanghaiforScienceandTechnology,Shanghai200093,China; 2.ShanghaiKeyLaboratoryofMultiphaseFlowandHeatTransferinPowerEngineering,UniversityofShanghaiforScienceandTechnology,Shanghai200093,China)

Numerical simulations were performed based on the open source computational fluid dynamics (CFD) software OpenFOAM.The numerically obtained velocity profiles over the dune were compared with the field measured results to validate the accuracy of the numerical simulation.On this basis,the velocity profiles at different heights over the dune as well as the velocity vector field were obtained.The sensitivity of the numerical simulation result to the computational mesh was discussed.Furthermore,the effect of the aspect ratio of the dune on the flow behavior was examined.The results show that the standardk-εturbulence model in OpenFOAM can be used to predict the velocity shift behavior of the flow field over the low sinusoidal dune.In the inner layer close to the dune surface,a phase advance of velocity curve relative to the dune profile is observable within a certain distance over the dune.As the height increases,the phase advance decreases and the velocity begins to flatten.The mesh number has little effect on the velocity profile in the outer layer,while it affects greatly on the velocity profile in the inner layer.It is worth noting that sufficient mesh number is required to attain accuracy information of the flow field over the dune. When the dune width keeps constant,with the decrease in the width-to-height ratio,the phase advance of the velocity profile over the dune increases and the recirculation zone could appear at the leeward.It is found that when the width-to-height ratio becomes smaller,the length and the intensity of the recirculation zone are greater.

computationalfluiddynamics(CFD);sinusoidaldune;turbulentflow;velocityshift;meshsensitivity

1007-6735(2017)04-0313-07

10.13255/j.cnki.jusst.2017.04.002

2016-11-24

國家自然科學(xué)基金資助項(xiàng)目(51206113,51176128,51576130);上海市科委科研計(jì)劃項(xiàng)目(13DZ2260900)

王 碩(1992-),男,碩士研究生.研究方向:多相流數(shù)值模擬.E-mail:598509460@qq.com

凡鳳仙(1982-),女,副教授.研究方向:多相流動(dòng)與傳熱.E-mail:fanfengxian@hotmail.com

O 355

A