祁紫薇，王 碩，凡鳳仙，胡曉紅

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

沙丘是不同粒徑的顆粒在風(fēng)力的長期作用下堆積而形成的自然地形地貌的一種，其形態(tài)、大小各不相同。沙丘的形成和演變過程強(qiáng)烈受制于風(fēng)沙兩相流運(yùn)動(dòng)規(guī)律，與氣流特征、顆粒運(yùn)動(dòng)緊密聯(lián)系，是一種極其復(fù)雜的物理過程。在研究沙丘形成和演變動(dòng)力學(xué)時(shí)應(yīng)充分考慮沙丘外形、沙丘表面流場、沙粒粒徑等多方面因素。其中，沙丘表面流場決定著局地沙粒輸運(yùn)的強(qiáng)度，在沙丘演變動(dòng)力學(xué)研究中占據(jù)重要地位。在沙丘表面流場的3 類研究方法中，野外觀測研究對象單一，而且受到惡劣的觀測環(huán)境以及儀器精度的限制，加之背風(fēng)坡的湍流特性十分復(fù)雜，依靠野外觀測難以對背風(fēng)坡流場結(jié)構(gòu)進(jìn)行精確描述；風(fēng)洞實(shí)驗(yàn)存在同時(shí)保證幾何相似和動(dòng)力學(xué)相似的困難；數(shù)值模擬能夠方便、靈活地調(diào)整計(jì)算參數(shù)，便于獲得流場結(jié)構(gòu)的細(xì)節(jié)信息，是研究沙丘表面流場特性的一種有效路徑[1-2]。

在沙丘表面流場數(shù)值模擬研究中，計(jì)算流體力學(xué)軟件PHOENICS[3]、FLUENT[4-6]、OpenFOAM[2,5,7-8]等以及自編程序[9]先后得到應(yīng)用，所采用的描述流場的數(shù)學(xué)模型包括標(biāo)準(zhǔn)k-ε、RNGk-ε、k-ωSST 湍流模型、大渦模擬等，文獻(xiàn)[2]對此作了較為詳盡的評述。由于不同的研究者所針對的沙丘形貌和尺度差別很大，采用的計(jì)算軟件及湍流模型也有所不同，所得的結(jié)論缺乏普適性。此外，已有數(shù)值模擬研究往往基于野外觀測或風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)考量模型的準(zhǔn)確性，對沙丘表面流場特性，特別是正弦形沙丘背風(fēng)坡回流區(qū)特性的研究仍不夠。筆者所在課題組曾借助OpenFOAM 的標(biāo)準(zhǔn)k-ε模型對正弦形沙丘表面流場開展數(shù)值模擬，將模擬結(jié)果與實(shí)驗(yàn)結(jié)果相結(jié)合驗(yàn)證了數(shù)值模擬的正確性和網(wǎng)格敏感性，給出了沙丘表面內(nèi)層和外層的速度分布曲線。在此基礎(chǔ)上，本文采用與文獻(xiàn)[2]一致的k-ε模型，通過改變沙丘的高度和寬度，開展數(shù)值模擬，考察沙丘表面流場的結(jié)構(gòu)與背風(fēng)坡回流區(qū)特性，以確定正弦形沙丘幾何參數(shù)對其表面背風(fēng)坡回流區(qū)的定量影響，為進(jìn)一步研究沙丘演變動(dòng)力學(xué)提供基礎(chǔ)和參考。

1 計(jì)算模型

描述湍流的模型，如標(biāo)準(zhǔn)k-ε、RNGk-ε、k-ωSST 湍流模型、大渦模擬等各有特點(diǎn)。標(biāo)準(zhǔn)k-ε模型是工程流場計(jì)算中主要的工具，其適用范圍廣，具有合理的計(jì)算成本和較高的計(jì)算精度，在流場模擬中具有廣泛的應(yīng)用。RNGk-ε模型相比標(biāo)準(zhǔn)k-ε模型修正了耗散率方程，對于復(fù)雜的剪切流，具有大應(yīng)變率、旋渦、分離等的流動(dòng)問題，具有比標(biāo)準(zhǔn)k-ε模型更高的精度。k-ωSST 模型包含了來自ω方程的交叉擴(kuò)散，在湍流粘度中考慮了湍流剪應(yīng)力的傳播，使得其用于流場模擬具有更高的精度。大渦模擬介于直接數(shù)值模擬（DNS）與雷諾平均納維-斯托克斯（RANS）方法之間，大尺度漩渦直接模擬，小尺度漩渦用RANS 方法求解，在實(shí)際工程應(yīng)用中需要很大的計(jì)算代價(jià)。由于在前期研究中，利用標(biāo)準(zhǔn)k-ε模型取得了與實(shí)驗(yàn)吻合良好的結(jié)果[2]，而標(biāo)準(zhǔn)k-ε模型計(jì)算成本相對較低，因此本文模擬中采用標(biāo)準(zhǔn)k-ε模型。

沙丘表面流場可用笛卡爾坐標(biāo)系下的雷諾平均納維-斯托克斯（RANS）方程進(jìn)行描述，可寫為

式中：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-ε湍流模型給出了湍動(dòng)能k和耗散率ε的方程如下：

式中：σk，σε，C1ε，C2ε為模型常數(shù)，在大氣邊界層模擬中通常采用σk=1.0，σε=1.30，C1ε=1.21，C2ε=1.92[10-12]；μt為湍流粘度；Eij為平均應(yīng)變率張量。

式中，模型常數(shù)Cμ=0.03[10-12]。

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

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

本文的研究是針對文獻(xiàn)[13]的沙丘形貌進(jìn)行的，即沙丘為狹長型沙丘（長30 m），如果按實(shí)際尺寸進(jìn)行模擬，計(jì)算網(wǎng)格數(shù)目龐大，導(dǎo)致計(jì)算成本很高。由于沙丘表面速度變化主要集中在來流方向（x向）和豎直方向（z向），因而可忽略y向的速度分量，開展二維模擬。事實(shí)上，已有研究也表明二維模擬能夠得到和實(shí)驗(yàn)吻合良好的結(jié)果[2]。圖1 給出了二維計(jì)算區(qū)域示意圖，沙丘位于計(jì)算區(qū)域底面中心，其輪廓曲線為

式中：zd為沙丘表面的z坐標(biāo)；A為正弦曲線的振幅；k為波數(shù)，k=2π/λ，λ為波長。

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

數(shù)值模擬時(shí)，為避免邊界效應(yīng)，進(jìn)口、出口、頂面與沙丘的距離應(yīng)足夠大。為此，計(jì)算區(qū)域設(shè)置時(shí)取Δx=4.35λ，Δz=18.20A，研究中發(fā)現(xiàn)，如果增大Δx和Δz，計(jì)算結(jié)果幾乎不發(fā)生變化。

2.2 網(wǎng)格劃分

為減少計(jì)算量，在劃分網(wǎng)格時(shí)，可以對于速度變化較大的區(qū)域采用較密的網(wǎng)格，對于速度變化較小的區(qū)域采用較疏的網(wǎng)格?；诖?，將由Δx和Δz確定的沙丘表面流場計(jì)算區(qū)域劃分為3 個(gè)子區(qū)域，如圖2 所示。其中，區(qū)域III 的大小由=λ×（4.55A）確定，區(qū)域II 的大小由（1.09λ）×（9.10A）確定。確定網(wǎng)格大小的原則為：區(qū)域I 的網(wǎng)格最大，為0.4 m×0.4 m；區(qū)域Ⅱ的網(wǎng)格邊長為區(qū)域I 的1/2，區(qū)域Ⅲ的網(wǎng)格邊長為區(qū)域Ⅱ的1/2。此外，在計(jì)算區(qū)域底面（即包含沙丘的地面）添加5 層邊界層，在遠(yuǎn)離下表面的方向上，邊界層網(wǎng)格厚度成比例增加，比例系數(shù)為1.1。

圖2 網(wǎng)格劃分所采用的3 個(gè)區(qū)域Fig.2 Three regions used for mesh generation

2.3 邊界條件

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

式中：κ為馮·卡門常數(shù)，κ=0.41；u*為剪切速度；z0為空氣動(dòng)力學(xué)粗糙長度。數(shù)值模擬中，選擇u*=0.5 m/s，z0=0.83 mm。

進(jìn)口湍動(dòng)能k和耗散率ε的表達(dá)式為[14-15]

底面采用壁面邊界條件，頂面設(shè)置為移動(dòng)壁面，頂面風(fēng)速由式（8）確定，出口設(shè)置為壓力出口。此外，基于OpenFOAM 進(jìn)行二維數(shù)值模擬時(shí)，需要在y向采用單層網(wǎng)格，并將相應(yīng)的邊界設(shè)置為“empty”。

2.4 求解控制

數(shù)值模擬中，采用SIMPLE 算法[16-18]對壓力和動(dòng)量方程進(jìn)行耦合求解，收斂條件為質(zhì)量和動(dòng)量通量的殘差達(dá)到10-8，壓力的殘差達(dá)到10-7；同時(shí)，壓力的松弛因子設(shè)置為0.3，速度、湍動(dòng)能和耗散率的松弛因子均設(shè)置為0.7；所采用的插值格式見表1。

表1 數(shù)值模擬采用的插值格式Tab.1 Interpolation schemes used in numerical simulations

3 結(jié)果與討論

3.1 沙丘高度對流場的影響

正弦形沙丘的振幅反映了沙丘的高度，為探討沙丘高度對沙丘表面流場的影響，保持波長λ=46 m 不變，在振幅A=1.65，2.20，2.75，3.30 m情況下開展數(shù)值模擬。模擬得到的沙丘表面不同高度處的速度曲線如圖3 所示，其中圖（a）～（d）分別對應(yīng)于上述4 種振幅情況，u0為對應(yīng)曲線上最大速度和最小速度的平均值。由圖3（a）可以看出：在沙丘高度較小時(shí)，如A=1.65 時(shí)，在距沙丘表面較高的高度范圍內(nèi)（1.5～5.0 m），速度的分布與沙丘外形類似，即呈對稱形，速度最大值出現(xiàn)在沙丘坡頂位置；在貼近沙丘表面的區(qū)域（沙丘表面0.5 m 高度以內(nèi)區(qū)域），速度最大值出現(xiàn)在沙丘迎風(fēng)坡，即坡頂上游，并且受沙丘輪廓影響，背風(fēng)坡流速先降低后增加，但速度值始終為正值，這表明在沙丘的背風(fēng)坡無明顯回流區(qū)。對比圖3（a）～（d）可知，隨著沙丘高度的增加，在A=2.20 m 時(shí)，貼近沙丘表面的區(qū)域中背風(fēng)坡流速出現(xiàn)負(fù)值，表明背風(fēng)坡流場出現(xiàn)了回流區(qū)。圖中，流速為負(fù)值的沙丘寬度（x）和高度區(qū)間反映了回流區(qū)的范圍，負(fù)的流速在數(shù)值上反映了回流區(qū)的強(qiáng)度?？梢?，隨著沙丘高度的增加，背風(fēng)坡回流區(qū)范圍變大、強(qiáng)度增強(qiáng)。圖3 也反映出與沙丘迎風(fēng)坡相比，背風(fēng)坡流速的變化更為劇烈，因此背風(fēng)坡沙粒受到的剪切力作用也更大；在沙丘背風(fēng)坡出現(xiàn)回流區(qū)時(shí)，沙丘表面沙粒還將受到渦流影響，因此背風(fēng)坡沙粒較迎風(fēng)坡沙粒更易于發(fā)生遷移。與此同時(shí)，貼近沙丘表面的氣流經(jīng)過沙丘坡頂后，在背風(fēng)坡轉(zhuǎn)而向下，同時(shí)水平方向速度迅速降低并出現(xiàn)回流，可以推測迎風(fēng)坡加速氣流挾帶的沙粒在輸運(yùn)到背風(fēng)坡時(shí)，傾向于降落在回流區(qū)的起始位置。

圖3 沙丘高度對沙丘表面流場的影響Fig.3 Effect of dune height on the flow field over the dune surface

3.2 沙丘表面回流區(qū)尺度的確定

沙丘背風(fēng)坡回流區(qū)的存在對沙丘的移動(dòng)以及演變有著重要影響，在對沙丘表面流場進(jìn)行定量研究時(shí)，回流區(qū)尺度（x向尺度用回流區(qū)長度表示，z向尺度用回流區(qū)高度表示）的確定是一個(gè)關(guān)鍵的問題。為確定回流區(qū)尺度，對數(shù)值模擬結(jié)果進(jìn)行處理時(shí)，首先繪制沙丘表面流場的流線圖，根據(jù)圖中氣流的回流與分離來確定回流區(qū)尺度。圖4 給出了不同高度沙丘的表面流場回流區(qū)。從圖4 可以看出，隨著沙丘高度的增加，回流區(qū)的長度和高度以及位置變化明顯，圖中示出了回流區(qū)的長度。結(jié)果表明：在振幅增加過程中，回流區(qū)長度增大；沙丘高度較小時(shí)，在沙丘背風(fēng)坡的底部貼近沙丘表面處形成扁平的連續(xù)渦流；隨著沙丘高度的增加，回流區(qū)中心從貼近沙丘表面的位置逐漸上移，回流區(qū)范圍逐漸增大。圖中所示3 種情況下，回流區(qū)長度與沙丘高度（2A）之比依次為1.02，3.33，4.23，圖4（b）和（c）中回流區(qū)長度范圍與Wiggs[19]的野外觀測數(shù)據(jù)相符，即回流區(qū)長度通常為沙丘高度的1.6～5.4 倍。但是，需要說明的是，野外觀測是對高度為9.6 m、坡頂與坡腳距離為86 m、寬度為130 m 的新月形沙丘進(jìn)行的。此外，江麗娟[6]對高度3～15 m，背風(fēng)坡坡度為32°的新月形沙丘表面流場進(jìn)行數(shù)值模擬研究，也得到回流區(qū)長度為沙丘高度的2～8 倍。本文針對的沙丘形狀與這些研究不同，但得出的回流區(qū)長度卻與之相符。

圖4 不同高度沙丘的回流區(qū)Fig.4 Recirculation zones of dunes with different heights

3.3 沙丘寬高比對回流區(qū)尺度的影響

圖5 給出了不同計(jì)算條件下回流區(qū)長度隨寬高比的變化關(guān)系。數(shù)值模擬時(shí)，采用了不同的計(jì)算區(qū)域設(shè)置，其中“原計(jì)算區(qū)域”表示按Δx=4.35λ、Δz=18.20A（見圖1）的原則確定的計(jì)算區(qū)域，“區(qū)域放大后”表示對計(jì)算區(qū)域長度和高度分別放大1.5 倍后的計(jì)算區(qū)域。對比相同沙丘幾何參數(shù)下“原計(jì)算區(qū)域”和“區(qū)域放大后”對應(yīng)的沙丘背風(fēng)坡回流區(qū)長度的計(jì)算結(jié)果，可知兩種計(jì)算區(qū)域下的計(jì)算結(jié)果吻合良好，顯示出本文數(shù)值模擬所選擇的計(jì)算區(qū)域?qū)τ谏城鸨砻媪鲌龌亓鲄^(qū)特性研究是合理的。沙丘輪廓線振幅A恒定時(shí)的曲線代表了回流區(qū)長度隨沙丘寬度的變化情況，可見：隨著沙丘寬度的增加，回流區(qū)長度減??；并且沙丘高度越大，回流區(qū)長度的減小越迅速。沙丘輪廓線波長λ恒定時(shí)的曲線代表了回流區(qū)長度隨沙丘高度的變化情況，可見，隨著沙丘高度的減小，回流區(qū)長度迅速減小?？偟卣f來，隨著寬高比的增加，回流區(qū)長度減小。然而，回流區(qū)長度不僅取決于寬高比，還取決于沙丘寬度與高度的數(shù)值。圖6 給出了回流區(qū)高度隨沙丘寬高比的變化情況?？梢?，沙丘寬高比越大，回流區(qū)高度越小。這表明在沙丘寬高比發(fā)生變化時(shí)，回流區(qū)在橫向和縱向的尺度同時(shí)發(fā)生變化，沙丘背風(fēng)坡回流區(qū)對沙丘幾何參數(shù)很敏感。

圖5 回流區(qū)長度隨沙丘寬高比的變化關(guān)系Fig.5 Relationship between the length of recirculation zone and the aspect ratio of the dune

圖6 回流區(qū)高度隨沙丘寬高比的變化關(guān)系Fig.6 Relationship between the height of recirculation zone and the aspect ratio of the dune

3.4 回流區(qū)成因分析

圖7 給出了λ=23 m、A=1.1 m 情況下，沙丘背風(fēng)坡的渦流分布情況。由于在沙丘背風(fēng)坡氣流受到的剪切應(yīng)力不均勻，造成氣流流動(dòng)不穩(wěn)定，容易形成渦流。沙丘背風(fēng)坡的回流區(qū)剛出現(xiàn)時(shí)極不穩(wěn)定，由一系列附著于沙丘表面的小渦組成，如圖7 所示。對于理想的平坦表面，定常流動(dòng)條件下，氣流為邊界層流動(dòng)，不發(fā)生流動(dòng)分離和回流。沙丘的存在，改變了氣流的流動(dòng)特性。對于正弦形沙丘，沙丘高度增加或?qū)挾葴p小造成沙丘坡度增大，對風(fēng)的阻礙作用增強(qiáng)，使氣流通過沙丘坡頂和背風(fēng)坡坡腳時(shí)，產(chǎn)生更大的壓力差，從而在背風(fēng)坡更易產(chǎn)生漩渦，形成回流區(qū)，進(jìn)而影響沙粒的遷移。需要說明的是，本文模擬著重從流場的角度探討了沙丘表面回流區(qū)特性，同時(shí)考慮沙粒的運(yùn)動(dòng)，利用氣固兩相流理論研究沙丘的演變過程將是下一步需要重點(diǎn)開展的工作。

圖7 沙丘（λ=23 m，A=1.1 m）背風(fēng)坡的渦流分布情況Fig.7 Eddy distribution at leeward side of the dune（λ=23 m,A=1.1 m）

4 結(jié)論

基于OpenFOAM 的標(biāo)準(zhǔn)k-ε模型，利用數(shù)值模擬方法研究了不同幾何參數(shù)條件下正弦形沙丘表面流場特性，特別是背風(fēng)坡回流區(qū)特性，給出了不同條件下回流區(qū)的特點(diǎn)，以及回流區(qū)長度和高度隨沙丘寬高比的變化關(guān)系，分析了回流區(qū)的成因，得到以下結(jié)論：

a.沙丘高度較低時(shí)，在背風(fēng)坡貼近沙丘表面的區(qū)域氣流速度先降低后增加；隨著沙丘高度的增加，背風(fēng)坡氣流速度降低更為迅速，氣流不穩(wěn)定，出現(xiàn)回流區(qū)。

b.在沙丘寬度或高度一定的情況下，隨著寬高比的增加，回流區(qū)尺度減小；但是回流區(qū)尺度不僅受沙丘寬高比的影響，還取決于沙丘寬度與高度的數(shù)值。

c.沙丘背風(fēng)坡回流區(qū)萌生時(shí)貼近沙丘表面出現(xiàn)一系列小渦，且分布不規(guī)律，隨著沙丘高度的增加或?qū)挾鹊臏p小，背風(fēng)坡氣流不穩(wěn)定性加劇，渦的尺度增大，強(qiáng)度增加，演變?yōu)榛亓鲄^(qū)。