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

(1.上海理工大學 能源與動力工程學院,上海 200093;2.上海理工大學 上海市動力工程多相流動與傳熱重點實驗室,上海 200093)

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

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

(1.上海理工大學 能源與動力工程學院,上海 200093;2.上海理工大學 上海市動力工程多相流動與傳熱重點實驗室,上海 200093)

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

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

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

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

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

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

2 計算模型

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

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

為使方程封閉,引入標準k-ε湍流模型計算湍動能k和耗散率ε,即

其中,

(5)

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

(6)

3 計算條件及方法

3.1 計算區(qū)域

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

zd=zref+Acos(Kx)

(7)

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

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

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

3.2 網(wǎng)格劃分

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

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

3.3 邊界條件

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

(8)

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

湍動能和耗散率的表達式為[20-21]

(9)

(10)

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

3.4 求解控制

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

4 結果與討論

4.1 數(shù)值模擬結果的驗證

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

4.2 沙丘表面流場速度分布

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

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

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

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

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

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

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

4.4 沙丘寬高比對流場的影響

為研究正弦形沙丘幾何參數(shù)對流場特性的影響,對不同寬高比的沙丘表面流場進行數(shù)值預測.在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處,隨著沙丘寬高比的減小,一方面,速度分布曲線相對于沙丘輪廓線的相位提前量增大,迎風坡氣流速度變化加劇;另一方面,當沙丘寬高比λ/(2A)≤10.5時,背風坡氣流速度出現(xiàn)負值,表明氣流流過坡頂后下沉形成回流區(qū).沙丘迎風坡的氣流加速和背風坡的回流區(qū)與沙丘侵蝕率有重要關聯(lián),沙丘寬高比越小,侵蝕越顯著,沙丘移動和演變越迅速.在距離沙丘表面0.5 m處,速度分布曲線與沙丘輪廓曲線的相位差均較小,尤其是在沙丘寬高比λ/(2A)≥10.5時,相位差幾乎為0;同時,沙丘寬高比越小,迎風坡氣流加速越快.沙丘寬高比λ/(2A)≥10.5時,在背風坡氣流速度全為正;而λ/(2A)≤8.36時,可觀察到速度負值,且高寬比越小,速度大小越偏離0.這些結果表明,隨著沙丘寬高比的減小,回流區(qū)長度和強度均增大.不同尺度沙丘表面流場特性研究是后續(xù)利用氣固兩相流數(shù)值模擬方法考察沙粒輸運,進而認識和掌握沙丘的形成機制、發(fā)展規(guī)律及變化趨勢的基礎.

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

5 結 論

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

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

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

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

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

[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] 呂萍,董治寶.沙丘動力學數(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].地理科學,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] 賈攀,王元,張洋.縮尺沙丘演化特性的風洞實驗研究[J].西安交通大學學報,2013,47(5):67-71.

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

[10] 周曉斯,王元,李志強.小尺度新月形沙丘背風側流場特性的大渦模擬分析[J].西安交通大學學報,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ù)值傳熱學[M].2版.西安:西安交通大學出版社,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

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

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

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

O 355

A