歐拉顆粒流模型下的近床面風沙躍移速度分布特征

劉 芳, 金阿芳, 熱依汗古麗·木沙

(新疆大學機械工程學院, 烏魯木齊 830047)

風沙顆粒在沙床表面的運動是一類復雜的多相流過程,隨著風速的增加,粒徑小于500 μm的沙粒最先被風吹起。當氣流對沙層表面施加的剪應力超過閾值0.05 N/m2時,沙粒就會發生躍移運動,相較于沙粒的蠕移和懸移運動,躍移運動所需風速最小,且是風沙輸移中沙量的主要來源,因此沙粒的躍移運動是風沙流研究中最重要的物理運動狀態。

從風沙邊界層動力學的研究角度來看,在大氣邊界層中,風沙流是固體顆粒與氣流的相互運動,能在大氣邊界層中形成明顯受風沙運動干擾的反饋兩相耦合作用的貼地氣流層。在風沙流中,由于受到氣相裹挾和湍動、沙粒運動和碰撞等因素的影響,風沙流動行為具有高度的隨機性、無規律性和結構不穩定性,即風沙流在躍移階段的運動具備稠密顆粒氣固兩相耦合的所有非線性系統特征。鄒學勇等[1]利用隨機概率方法推導出起躍沙粒垂直初速度的分布函數,進一步解釋了風沙流結構呈對數分布。亢力強等[2]采用相位多普勒粒子分析儀(phase Doppler particle analyzer，PDPA) 測量了風沙兩相流動中沙床面上沙粒起跳和碰撞速度的概率分布以及不同高度處沙粒速度的概率分布,實驗結果表明沙床面上沙粒起跳和碰撞速度概率分布均呈現對數正態分布規律, 沙粒垂直速度概率分布在不同高度處都可用正態分布函數來描述,碰撞和起跳角度均可表示為指數分布函數。肖鋒軍等[3]同樣使用PDPA測量了近床面沙粒的躍移速度,并給出了沙粒沖擊速度和角度的聯合概率分布函數。程旭[4]在對沙粒的法向速度、水平速度和起動角度的分布研究發現,各曲線均表現出正態分布的特征,且沙粒起動速度的頻率分布與風速無關。何麗紅等[5]對已有的沙粒起跳垂直速度分布函數進行了非參數檢驗,最后提出了躍移沙粒起跳垂直速度分布函數服從Weibull分布的結論。羅生虎等[6]對風沙流中沙粒的起跳初速度進行分析發現該速度分布函數服從瑞利分布。金阿芳等[7]運用光滑粒子流體動力學法(smooth particle hydrodynamics，SPH)方法模擬了平坦沙床上沙粒的起跳至回落過程,并對隨機抽取的1 000個沙粒進行了統計分析后發現,躍移沙粒起跳速度和起跳角的概率密度函數分別服從對數正態分布和指數分布。

為了完整考慮顆粒相的湍流輸運過程,并考慮現實環境中計算機硬件及CPU速度的限制,運用歐拉-歐拉兩相顆粒流模型對風沙流運動進行模擬。由于歐拉模型涉及的參數眾多,眼下有關歐拉-歐拉模型對風沙流的模擬研究還不充分,張偉等[8]利用歐拉雙流體模型研究了該模型中風速和碰撞恢復系數對風沙運動的影響。在此基礎上,運用Fluent軟件建立歐拉兩相流中的顆粒流模型,對風沙流的起動、發展及衰退過程進行模擬。沙粒的速度分布作為研究輸沙量分布的重要參數,主要包含水平和鉛垂兩個方向的運動信息。對風沙流整個運動流程中的物理數據進行統計分析,從沙相速度分布特征中研究沙粒躍移運動中的普遍規律。

1 風沙流數學模型

1.1 近地表層風沙運動的控制方程

在大氣邊界層里,風沙運動最活躍的區域是下部低于2 m的近地表層,受地表條件的影響,貼地層的速度梯度大并伴有明顯的湍流特性,流場受下墊面影響大。風沙活動取決于地表沙粒和大氣之間的動力和熱力作用,在非極端情況的天氣條件下,熱量因素與動力因素相比影響較小,可忽略不計。用歐拉-歐拉模型模擬風沙流運動時,氣體和沙粒相具有各自的質量和動量平衡輸運方程。各相質量平衡方程為

(1)

式(1)中:αk為體積分數；ρk為密度；Uk為k相的平均速度，下標k為g表示氣體,k為p表示顆粒相。

各相的動量平衡方程為

(2)

式(2)中:Pg為平均氣體壓力;gi為重力加速度;∑k,ij為有效應力張量;Ik,i為不考慮平均氣體壓力影響的相間平均動量傳遞速率。當只考慮相之間的拖曳力和平均粒子弛豫時間標度τp時,Ik,i為

(3)

式(3)中:

(4)

(5)

(6)

式中：CD為阻力系數,由于沙粒起動不僅受氣流平均速度的影響,湍流效應也會影響沙粒的擴散,引入顆粒雷諾數Rep來說明沙粒對氣流的影響；dp為顆粒直徑；vg為氣體分子運動黏度；局部瞬時相對速度vr,i等于局部粒子速度up,i與瞬時氣體速度ug,i之差；Vr,i為平均顆粒流相對速度,Vr,i=Up,i-Ug,i。顆粒應力張量∑p,ij定義為

(7)

式(7)中:碰撞壓力Pp和體積黏度λp決定了粒子間碰撞時的能量損失。

1.2 兩相的體積分數方程

多相流各相間的相互貫穿和融合用相體積分數的概念來描述,這里沙相和氣相體積分數分別表示為αp和αg,體積分數能表示各相所占的空間大小,其中每一相都分別滿足質量守恒定律和動量守恒定律。多相流中k相的體積Vk定義為

(8)

1.3 湍流模型

由于氣流表觀雷諾數的范圍為15 000～50 000,且底部沙粒形成的粗糙表面會導致氣流速度呈明顯的湍流分布。這進一步說明了在近地面高度的沙面流場屬于湍流邊界層。現運用標準κ-ε湍流模型對方程進行封閉。

2 數值計算

計算區域的設置如圖1所示。根據風沙流躍移運動中的參數變化特點,將計算區域設置為1 m×1 m的二維矩形區域,網格為均勻四邊形網格,數量為111 556個。沙粒密度ρp=2 650 kg/m3,空氣密度ρg=1.255 kg/m3,沙床厚度為0.05 m。

圖1 數值計算平坦沙床示意圖Fig.1 Schematic diagram of numerical calculation of flat sand bed

2.1 歐拉顆粒流模型設置

在歐拉模型中,將沙相設置為顆粒相,顆粒黏度選用syamlal-obrien模型,顆粒體積黏度選用lun-et-al模型,最大堆積率為0.63,曳力模型設置為gidaspow。

2.2 沙相設置

在矩形計算域中標記出0.05 m×1 m的沙相區域,沙相區域網格數為5 678個,流場初始化后利用Patch功能將沙相區域的體積分數設為1。

2.3 邊界條件

算例中邊界條件的設置為左邊界為空氣相凈風速度入口,下邊界無滑移壁面邊界,為使風場發展充分,且符合自然風條件,模型上邊界設置為對稱邊界,右邊界設置為壓力出口邊界。

2.4 速度入口邊界條件

3 計算結果

在粒徑為0.000 1 m球形沙粒組成的平坦靜止沙床上進行模擬,摩阻風速為0.265 m/s,該風速下,不同時刻在沙床中心沿高度y分布的沙相體積分數αp，如圖2所示。由圖2可見,除沙床及其表面附近區域外,從沙粒起動到穩定的不同時刻αp在lnαp-h的半對數坐標系中均呈直線分布,且αp在10-3～10-7范圍內的數據與實驗結果[9-10]吻合良好,這說明該模型能有效模擬風沙流的躍移運動,且相比于實驗結果其數據范圍更為廣泛,即該模型有一定的預測能力。在0.3 m高度處(距初始床面0.25 m),沙相濃度出現飽和,其中飽和高度的αp=10-7,隨著時間的推移,沙相在氣流中的濃度逐漸保持穩定,而在飽和層上方,沙相體積分數衰減加快。不同時刻沙相體積分數云圖如圖3所示。

圖2 不同時刻沙相體積分數隨高度變化Fig.2 The volume fraction of sand phase changes with height at different times

圖3 不同時刻沙相體積分數云圖Fig.3 Cloud maps of sand phase volume fraction at different times

沙粒的速度分布是研究輸沙量分布的重要參數,與拉格朗日法不同,歐拉方法不能得到單顆沙粒的運動速度,但通過對網格內沙相速度的平均就可以得到不同時刻沙相的速度信息,在此涉及的沙粒速度均為網格節點處沙相的平均速度。在對全流程不同高度處沙粒垂向速度進行統計發現,沙粒垂向速度在不同高度處均呈正態分布,該結論與文獻[2]結果相一致。由圖4可以看出, 沙粒垂向速度分布在不同高度均存在高度集中區域,且中心隨高度向右偏移,隨著高度的增加,沙粒垂向速度分布越發松散,這是由于沙粒濃度及沙粒與氣流速度兩相發生耦合引起的反饋。

圖4 全流程不同高度處沙粒垂向速度分布Fig.4 Vertical velocity distribution of sand grains at different heights in the whole process

現有研究的沙粒速度分布是通過統計若干顆沙粒的速度信息后得到的,統計結果受沙粒數目影響大,且常用的實驗方法在統計沙粒垂向速度時,往往會遺漏掉速度為0的沙粒,這樣得到的沙粒速度分布是不準確的。根據不同時刻沙粒體積分數αp的變化,對αp≥10-7的沙粒速度信息進行統計,得到不同時刻沙床附近所有沙粒在躍移運動中的速度分布,如圖5所示。

由圖5可知，不同時刻αp≥10-7區域的沙相水平速度呈對數正態分布。進一步分析得到,沙相水平速度在沿流場方向均是正值,且速度分布除0時刻外整體向左偏移,分布并不對稱,因此沙相水平速度并不符合普通正態分布,而用對數正態分布擬合較好。即lnvx～N(μ,σ2),其中沙相水平速度vx為取值為正的連續隨機變量,vx的概率密度為

f(vx,μ,σ)=

(9)

式(9)中:μ、σ2分別為原始數據的均值和方差,它們的值與風沙顆粒的位置、形狀及尺度參數相關。由圖5可知,在不同時刻沙粒的水平速度分布并不相同,隨時間推移,vx的均值變小,且整體分布趨于平緩。

圖5 不同時刻沙相水平速度分布Fig.5 Horizontal velocity distribution of sand phase at different times

圖6所示為在不同時刻,αp≥10-7區域的沙相垂向速度vy分布。圖6中的分布曲線為拉普拉斯分布(雙指數分布)。可以看出,沙粒的上升段(vy>0)和下降段(vy<0)速度都服從指數分布,在沙粒初始運動階段,起跳沙粒數量逐漸增多,隨后沙粒開始降落,參與運動的沙粒數目逐漸增大,最終保持穩定的跳躍狀態,但沙粒的上升和下降速度分布并不對稱(分布向下降段偏移),盡管沙粒對床面的沖擊和碰撞會激發更多沙粒參與運動,但沙粒對風速的阻礙會導致風速減小,若沒有風力和沙源的補充,可以預測沙面最終能趨于平靜,即再無沙粒運動。沙相垂向速度分布可以表示為vy～La(μ,b),vy的概率密度為

圖6 不同時刻沙相垂向速度分布Fig.6 Vertical velocity distribution of sand phase at different times

(10)

式(10)中：μ為位置參數；b>0為尺度參數。

4 結論

采用二維歐拉顆粒兩相流模型對風沙的流動過程進行模擬,實現了風沙流沙相體積分數及躍移層沙粒速度分布的研究,得到以下結論。

(1)風沙流從起動到穩定的過程中沙相體積分數在lnαp-h的半對數坐標系中均呈直線分布,該模型不僅能有效模擬風沙流的躍移運動,還具有一定的預測能力。沙相濃度在0.3 m高度處出現飽和,其中飽和高度的沙相體積分數為10-7,在飽和層上方,沙相體積分數衰減加快。

(2)通過對體積分數大于等于10-7的沙粒速度信息進行統計,得到不同時刻沙相水平速度呈對數正態分布、沙相垂向速度分布為拉普拉斯分布(雙指數分布),并給出了相應的速度分布概率密度函數。

(3)盡管風沙流是一種多耦合、非線性的隨機運動,但利用歐拉顆粒流模型及數理統計的方法能夠有效地揭示風沙顆粒流體的動力學運動特征。由于沒有考慮混合粒徑及沙粒形狀對速度分布的影響,這些因素會影響概率密度函數中參數的確定,這也將是后期研究的一個方向。