正弦形沙丘背風坡回流區特性研究

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

（1.上海理工大學 能源與動力工程學院，上海 200093；2.上海理工大學 上海市動力工程多相流動與傳熱重點實驗室，上海 200093）

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

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

1 計算模型

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

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

式中：t為時間；u為流速；p為壓力；ρ為空氣密度；μ為空氣動力粘度；δij為克羅內克函數，當i=j時，δij=1，當i≠j時，δij=0。

為使方程封閉，標準k-ε湍流模型給出了湍動能k和耗散率ε的方程如下：

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

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

2 計算條件及方法

2.1 計算區域

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

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

圖1 計算區域示意圖Fig.1 Schematic diagram of the simulation domain

數值模擬時，為避免邊界效應，進口、出口、頂面與沙丘的距離應足夠大。為此，計算區域設置時取Δx=4.35λ，Δz=18.20A，研究中發現，如果增大Δx和Δz，計算結果幾乎不發生變化。

2.2 網格劃分

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

圖2 網格劃分所采用的3 個區域Fig.2 Three regions used for mesh generation

2.3 邊界條件

進口采用速度邊界條件，速度值可由大氣邊界層風速的對數分布率給出[13]

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

進口湍動能k和耗散率ε的表達式為[14-15]

底面采用壁面邊界條件，頂面設置為移動壁面，頂面風速由式（8）確定，出口設置為壓力出口。此外，基于OpenFOAM 進行二維數值模擬時，需要在y向采用單層網格，并將相應的邊界設置為“empty”。

2.4 求解控制

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

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

3 結果與討論

3.1 沙丘高度對流場的影響

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

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

3.2 沙丘表面回流區尺度的確定

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

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

3.3 沙丘寬高比對回流區尺度的影響

圖5 給出了不同計算條件下回流區長度隨寬高比的變化關系。數值模擬時，采用了不同的計算區域設置，其中“原計算區域”表示按Δx=4.35λ、Δz=18.20A（見圖1）的原則確定的計算區域，“區域放大后”表示對計算區域長度和高度分別放大1.5 倍后的計算區域。對比相同沙丘幾何參數下“原計算區域”和“區域放大后”對應的沙丘背風坡回流區長度的計算結果，可知兩種計算區域下的計算結果吻合良好，顯示出本文數值模擬所選擇的計算區域對于沙丘表面流場回流區特性研究是合理的。沙丘輪廓線振幅A恒定時的曲線代表了回流區長度隨沙丘寬度的變化情況，可見：隨著沙丘寬度的增加，回流區長度減小；并且沙丘高度越大，回流區長度的減小越迅速。沙丘輪廓線波長λ恒定時的曲線代表了回流區長度隨沙丘高度的變化情況，可見，隨著沙丘高度的減小，回流區長度迅速減小。總地說來，隨著寬高比的增加，回流區長度減小。然而，回流區長度不僅取決于寬高比，還取決于沙丘寬度與高度的數值。圖6 給出了回流區高度隨沙丘寬高比的變化情況。可見，沙丘寬高比越大，回流區高度越小。這表明在沙丘寬高比發生變化時，回流區在橫向和縱向的尺度同時發生變化，沙丘背風坡回流區對沙丘幾何參數很敏感。

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

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

3.4 回流區成因分析

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

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

4 結論

基于OpenFOAM 的標準k-ε模型，利用數值模擬方法研究了不同幾何參數條件下正弦形沙丘表面流場特性，特別是背風坡回流區特性，給出了不同條件下回流區的特點，以及回流區長度和高度隨沙丘寬高比的變化關系，分析了回流區的成因，得到以下結論：

a.沙丘高度較低時，在背風坡貼近沙丘表面的區域氣流速度先降低后增加；隨著沙丘高度的增加，背風坡氣流速度降低更為迅速，氣流不穩定，出現回流區。

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

c.沙丘背風坡回流區萌生時貼近沙丘表面出現一系列小渦，且分布不規律，隨著沙丘高度的增加或寬度的減小，背風坡氣流不穩定性加劇，渦的尺度增大，強度增加，演變為回流區。