近床面風沙流的顆粒擬流體大渦模擬分析

周曉斯，王元，李志強

（西安交通大學流體機械及工程系，710049，西安）

目前，針對風沙流模擬的研究大多集中于以對沙相采用拉格朗日坐標、對氣相采用歐拉坐標而建立的相關(guān)離散顆粒模型進行的數(shù)值計算［1-2］?；陔x散顆粒模型的風沙流數(shù)值計算，因其可以獨立追蹤單顆沙粒運動，所以在研究風沙相間作用及沙粒間碰撞作用等微觀現(xiàn)象方面具有獨特的優(yōu)勢［2］。當顆粒濃度較大時，強烈的相間耦合使得流場的收斂變得困難［3］。計算表明，離散顆粒模型計算耗時要比擬流體模型高出4個數(shù)量級［4］。

劉大有等指出，顆粒運動研究只能采用連續(xù)介質(zhì)力學方法［5］。近床面大量沙粒躍移表現(xiàn)出了群體運動特性，通過開展基于連續(xù)介質(zhì)假設(shè)的擬流體模擬，能夠獲得類似于流體相的顆粒群宏觀參量。基于此，建立針對氣沙兩相擬流體模型來模擬風沙具有一定的理論意義和工程應(yīng)用價值。基于擬流體模型的兩相流模擬，在氣固本構(gòu)方程的建立及反映氣固相互作用的曵力模型適用性確定方面，目前并沒有形成一個普適的模型理論體系［6］。針對風沙運動的擬流體數(shù)值計算，一些研究者做出了有益嘗試［7-9］。然而，沙粒相在風場中體現(xiàn)出顯著的三維運動特征［2］，簡化的二維計算模型必然與真實運動特性存在較大偏差，同時將體現(xiàn)沙粒湍流脈動特性的顆粒擬溫度簡化為常數(shù)［9］，其計算結(jié)果的準確性尚需相關(guān)實驗驗證。

本文基于氣相湍流大渦數(shù)值模擬，將風沙流中沙粒相處理為擬流體相，引入顆粒動理學理論建立顆粒相本構(gòu)方程，采用自編程序?qū)炒裁孳S移層內(nèi)風沙流開展三維數(shù)值模擬，從而獲得了近床面沙粒群躍移運動的速度場、濃度場及擬溫度場等相關(guān)參數(shù)的空間宏觀分布特性；基于流場參數(shù)的統(tǒng)計特性，分析了躍移層內(nèi)反映沙粒湍流脈動特性的沙相脈動均方根速度參數(shù)變化及碰撞彈性恢復(fù)系數(shù)；基于模擬結(jié)果的對比分析，對所建模型的適用性及改進性進行了討論。

1 數(shù)學模型

1.1 控制方程

借鑒文獻［4］，改進后適用于本文模擬的單固體相大渦模擬氣沙兩相控制方程如下。

連續(xù)性方程

式中：k表示相，氣相時k為g，沙相時k為s；αk為相體積分數(shù)。

氣相動量方程

式中：μgt為氣相亞格子渦黏系數(shù)。本文基于大渦模擬Smagorinsky SGS模型，則有

式中：Δ 為濾波尺度，Δ=（Δx×Δy×Δz）1／3；D 為引入修正函數(shù)壁面模型的近壁衰減函數(shù)［10］，D=1-為Smagorinsky常數(shù)，計算時Cg=0.10。

沙相動量方程

式中：ps為濾波后的沙相壓強；μs為沙相剪切黏度；ξs為沙相表觀黏度；e為顆粒碰撞彈性恢復(fù)系數(shù)，實際沙粒是以非彈性碰撞形式耗散能量的，計算中e=0.70［1］；θ為顆粒擬溫度反映顆粒相速度脈動強弱，不同于顆粒本身的熱力學溫度；g0為徑向分布函數(shù)

其中αs，max為顆粒自由堆積時最大體積分數(shù)，針對沙相屬性，計算中αs，max=0.75。

針對顆粒擬溫度，建立如下方程

式中：Ks為顆粒熱導(dǎo)率，反映顆粒碰撞引起的能量交換

γs為顆粒脈動能量引起的耗散率

Φs為氣固兩相之間的脈動能量交換率

式（2）、（6）、（14）中β為相間曵力系數(shù)，是耦合氣相和沙相之間動量交換的關(guān)鍵參數(shù)。本文嘗試采用下面2種模型進行對比分析。

（1）Gidaspow 修正模型［11］，即

式中：轉(zhuǎn)換函數(shù)

（2）Felice曵力模型［12］，即

以上兩模型中顆粒雷諾數(shù)

1.2 數(shù)值方法

氣固兩相控制方程的求解采用同位網(wǎng)格下的SIMPLE算法，來處理兩相流中速度與壓強耦合的問題。兩相動量方程及顆粒擬溫度方程的對流項和擴散項的空間離散分別采用QUICK格式和中心差分格式；對于時間離散，對流項采用二階Adams-Bashforth格式，擴散項采用二階Crank-Nicolson格式；顆粒相體積分數(shù)方程采用QUICK空間格式及全隱時間格式。為了加速兩相速度的計算收斂，兩相動量方程采用加速收斂的PEA（Partial-Elimi-nation-Algorithm）技術(shù)［13］進行離散處理。數(shù)值計算流程如圖1所示。

圖1 數(shù)值計算流程圖

以上計算方法與文獻［13］采用的計算方法略有不同，區(qū)別在于本文采用顆粒相體積分數(shù)修正方程來修正最終的顆粒相速度。計算表明，采用此方法，計算迭代收斂速度及穩(wěn)定性均得到了顯著改善［14］。

2 計算模型及計算條件

2.1 計算模型

近床面風沙流場三維長方體計算域結(jié)構(gòu)及網(wǎng)格劃分如圖2所示。計算域在流向、垂向和展向上采用的空間尺度為Lx×Ly×Lz=0.10πm×0.10m×0.05πm，網(wǎng)格數(shù)為Nx×Ny×Nz=48×64×64，垂向沙床面進行網(wǎng)格加密，以提高計算精度。

2.2 計算條件

為成功模擬湍流流場隨時間的發(fā)展過程，風沙流的兩相大渦數(shù)值模擬中的計算域氣體相采用符合對數(shù)律分布的平均速度剖面，并疊加一定強度的速度脈動量對氣相速度場進行初始化。模擬流場環(huán)境溫度為常溫25℃，整個模擬過程中保持溫度不變。空氣密度ρg=1.225kg／m3，空氣動力黏度μgl=17.85μPa·s，沙相密度ρs=2 650kg／m3，沙相（簡化為球形顆粒）直徑ds=110μm。

圖2 計算域結(jié)構(gòu)簡圖及網(wǎng)格

計算域流向和展向均采用周期性邊界條件。在計算域法向上，沙床面處氣相速度采用無滑移條件，沙床面頂部采用對稱邊界條件。計算時初始壓力場為0。借鑒文獻［2］，計算時間步長為0.01ms。

沙床面沙相邊界條件采用Johnson和Jackson提出的壁面邊界條件［15］。床面相切的沙相速度

顆粒擬溫度

式中：ψ為鏡面反射系數(shù)，當沙粒與沙床面為完全彈性碰撞時ψ=1，為完全非彈性碰撞時ψ=0，由床面沙相屬性計算時ψ=0.01；ew為沙粒的壁面碰撞彈性恢復(fù)系數(shù)。近床面沙相濃度越高，ew越小，計算床面沙相屬性時ew=0.70。

式（18）中等號右端第二項為沙相速度滑移引起的產(chǎn)生項；式（19）為近床面沙粒非彈性碰撞引起的耗散項，與床面垂直的沙相速度分量均為0。

實質(zhì)的風沙流，在垂直于地表一定高度y以下沙粒群會因稠密堆積而不再運動。本文將不再運動的沙粒層設(shè)置為計算域下邊界靜止不動的壁面，即靜止沙床面，這種簡化設(shè)置有別于實際情況，但對研究躍移沙粒運動的本質(zhì)影響不大。t=0s時，在y＜10cm空間隨機設(shè)置一定的沙相濃度（沙子的體積分數(shù)）來表征初始時刻該空間位置處存在的沙粒數(shù)量，沙相初始速度為0。如果計算域宏觀的總沙相濃度隨時間不變，則認為流場達到動態(tài)穩(wěn)定態(tài)。

3 計算結(jié)果及分析

3.1 來流風場分析

圖3 凈風場時均風速廓線

3.2 沙相流場分析

特定工況下躍移層風沙流實驗表明，在沙相躍移底層（y＜2cm），對于ds=100～125μm的沙粒，氣沙兩相間的平均流向速度相對值Urel的變化范圍在1.162～1.660m／s之間［16］。為便于比較，取ds=110μm、Urel=1.50m／s計算2種曵力模型的曵力系數(shù)隨沙相濃度的變化，結(jié)果如圖4所示。由圖4可見，以沙相濃度作為唯一變量計算出的不同模型曵力系數(shù)均隨沙相濃度的增大而增大。在整個沙相濃度變化范圍內(nèi)，2種曵力模型預(yù)測出的曵力系數(shù)吻合度較高。0.20＜αs＜0.65時，Gidaspow修正模型的曵力系數(shù)略高于Felice曵力模型，在極濃（αs≈0.70）和極?。é羢≈0.00）情況下，2種曵力模型的曵力系數(shù)趨于一致。

圖4 曵力系數(shù)隨沙相濃度的變化

2種曵力模型下躍移層沙相流向速度的統(tǒng)計平均值沿y的變化如圖5所示。為便于對比分析，圖中還給出了Liu等風洞實驗測得的沙相平均流向速度統(tǒng)計值的變化［17］。由圖5可見，躍移層內(nèi)，2種曵力模型在近床面的沙相速度均比較小，隨著y的增加，速度不斷提高。這是躍移沙粒不斷被來流風場加速的緣故。在躍移中下層，沙相速度沿y變化的模擬值與實驗值較為接近，在躍移上層二者偏差較大，模擬值大于實驗值。文獻［17］的實驗數(shù)據(jù)表明，相同自由來流風速下，躍移沙相平均流向速度隨沙粒粒徑的增大而下降。在近床面，沙粒與床面以及沙粒之間的碰撞主導(dǎo)著沙粒運動；在離開床面較高處，來流風場成為沙粒運動的主導(dǎo)者，沙粒粒徑越小，沙粒與流體的跟隨性越強，沙粒易被加速，沙粒受風加速的時間越長，沙粒粒徑的影響導(dǎo)致的偏差凸顯。以上規(guī)律可以解釋本文模擬采用的ds=110μm單一沙粒粒徑下計算得到的在躍移層較高處的沙相平均流向速度，及Liu等在100μm＜ds≤200μm混合沙粒粒徑下得到的速度實驗值出現(xiàn)較大偏差的原因。由圖5還可見，2種曵力模型下計算得到的速度沿y的變化達到高度吻合。這是因為計算時風沙流中的沙相稀疏，該狀態(tài)下2種曵力模型的系數(shù)相當，所以2種模型的模擬結(jié)果達到一致。沙相平均流向速度沿y變化的模擬結(jié)果可用冪函數(shù)表示，即us=ayb，擬合回歸系數(shù)a=35.093 9，b=0.529 2，相關(guān)系數(shù)平方R2=0.995 71。

圖5 沙相平均流向速度沿y的變化

圖6 沙相濃度沿y的變化

2種曵力模型下躍移層沙相濃度統(tǒng)計平均值沿y的變化如圖6所示。由圖6可見，與Liu等實驗結(jié)果相比［17］，2種曵力模型的沙相濃度在躍移中下層的模擬值略高，在躍移上層偏低。這是因為本文模擬時在計算域流向和展向上均設(shè)定了周期性邊界條件，這2個方向上被夾帶流出區(qū)域的沙粒又由各個方向的入口返回，所以整個模擬過程中計算域平均沙相濃度保持不變。如果沙相濃度在躍移中下層被高估，則在躍移上層必然被低估。導(dǎo)致躍移上層出現(xiàn)偏差的另一種原因是，這2種曵力模型均低估了曵力系數(shù)，使得沙粒加速度偏小，沙粒弛豫時間加長，大量沙粒不能被來流風夾帶到高層，從而出現(xiàn)了躍移層上部沙相濃度的模擬值偏小的現(xiàn)象。沙相濃度沿y變化的模擬結(jié)果可用指數(shù)函數(shù)表示，即αs=aexp（by），擬合回歸系數(shù)a=9.427 4×10-5，b=-25.135 3，相關(guān)系數(shù)平方R2=0.999 97。

將本文Felice曵力模型下的沙相速度及濃度分布與Liu等的實驗進行對比可見，在相同來流風速及合理粒徑條件下，擬流體模擬得到的沙相速度及濃度分布在宏觀統(tǒng)計值上均與實驗值存在一定的偏差。為了實現(xiàn)氣沙兩相流場擬流體模型的準確模擬，需要對風沙流特性的擬流體氣固本構(gòu)封閉方程和相間作用模型進行改進。

定義沙相流向脈動均方根速度（RMS）［18］

圖7 沙相uRMS沿y的變化

2種曵力模型下躍移層內(nèi)uRMS沿y的變化如圖7所示。由圖7可見，躍移層內(nèi)模擬得到的uRMS的變化范圍在0.20～0.80m／s之間。風沙兩相流中，近床面的uRMS主要由剪切力產(chǎn)生，沿沙床面由高到低先漸增，后急劇減小，在近床面附近出現(xiàn)最大峰值。圖7顯示出了Kang等風洞實驗測得的uRMS沿y的變化［19］。實驗中混合沙粒粒徑變化范圍為170μm≤ds≤300μm，自由來流風速U0=9.3m／s。通過對比可知，模擬得到的uRMS沿y的分布規(guī)律與實驗結(jié)果基本一致，但模擬值均比實驗值小得多。這是模型模擬存在固有誤差、自由來流風速及沙粒粒徑不同的緣故。由圖7還可見，躍移頂層的模擬值有快速減小的趨勢。這可能是躍移頂層預(yù)測的沙相濃度偏?。ㄒ妶D6），沙粒之間碰撞作用顯著下降所致。為揭示沙粒粒徑對uRMS沿y的分布規(guī)律及模擬結(jié)果準確性的影響，在ds=110μm基礎(chǔ)上，基于實驗沙相粒徑分布范圍170μm≤ds≤300μm［19］，取粒徑ds分別為170、220、260、300μm 進行了計算。Felice曵力模型下ds對uRMS影響的計算結(jié)果如圖8所示。由圖8可見，不同ds下，近床面uRMS沿y均呈現(xiàn)出單一最大值，且隨沙粒粒徑的增大而增大。與實驗結(jié)果對比表明，實驗中uRMS沿y的分布介于模擬中最小ds和最大ds之間，實驗與模擬結(jié)果僅在躍移頂層偏差稍大。由于實驗采用混合粒徑，由此顯見，所建模型的模擬結(jié)果是正確、合理的。

圖8 Felice曵力模型下ds對uRMS影響的計算結(jié)果

擬流體模型中顆粒碰撞彈性恢復(fù)系數(shù)e是準確表征沙相運動特性的關(guān)鍵參數(shù)，顆粒之間碰撞的能量耗散取決于e值。本文計算中，對于給定的沙相屬性及速度變化范圍，e為常數(shù)。實質(zhì)上，e取決于顆粒之間碰撞的速度。Haff等指出，沙粒的外形無規(guī)則，所以準確的e很難通過實驗獲得，同時建議e的合理取值范圍在0.6～0.8之間［20］。本文計算、分析了e分別為0.6、0.7、0.8時躍移層內(nèi)uRMS沿y的變化。Felice曵力模型下e對uRMS的影響如圖9所示。由圖9可見：e分別為0.6、0.7、0.8時，uRMS在近床面對應(yīng)的單一峰值分別為0.702 9、0.799 3、0.961 3m／s；e較小，uRMS低，e較大，uRMS高。

圖9 Felice曵力模型下e對uRMS的影響

由圖9還可見，uRMS并非隨e的均勻增大呈現(xiàn)均勻遞增的趨勢，e越大，表明沙粒間碰撞的動能損失越小，攜帶較大動能余量的沙粒在風場中可以不斷地發(fā)生頻繁碰撞，最終導(dǎo)致沙相劇烈脈動。

4 結(jié) 論

（1）本文模型的模擬結(jié)果能較好地揭示近床面風沙流特性。躍移層內(nèi)，us隨y的增加呈遞增趨勢，躍移高層模擬結(jié)果與相關(guān)實驗存在較大偏差。

（2）躍移層內(nèi)，沙相濃度統(tǒng)計平均值隨y的增加呈遞減趨勢。與相關(guān)實驗結(jié)果對比可見，曵力模型低估了曵力系數(shù)，從而導(dǎo)致沙相濃度在躍移中下層的模擬結(jié)果偏高，而在躍移上層偏低。

（3）在近床面，躍移層內(nèi)的uRMS沿y分布呈現(xiàn)出典型的單一峰值，模擬結(jié)果的變化規(guī)律與實驗基本一致。在躍移高層，低估沙相濃度會導(dǎo)致uRMS隨y的增加快速減小。

（4）躍移層內(nèi)，uRMS隨沙相粒徑的增大而增大，模擬結(jié)果與實驗相符。uRMS隨e的增加而增大，e越大，沙粒間形成頻繁碰撞而導(dǎo)致沙相脈動越劇烈。

（5）在沙相濃度范圍內(nèi)，本文的2種曵力模型高度吻合。在躍移高層，模擬結(jié)果出現(xiàn)較大偏差，表明擬流體模型的優(yōu)勢主要體現(xiàn)在稠密氣固模擬上。通過改進氣相大渦模擬及沙相本構(gòu)方程，采用多粒徑混合沙粒的多相流風沙運動擬流體模擬，將是下一步開展的工作。

［1］ KANG L Q.Discrete particle model of aeolian sand transport：comparison of 2Dand 2.5Dsimulations［J］.Geomorphology，2012，139／140（15）：536-544.

［2］ 李志強，王元，王麗，等.風沙流中近床面沙粒三維運動的LES-DEM分析 ［J］.空氣動力學學報，2011，29（6）：784-788.LI Zhiqiang，WANG Yuan，WANG Li，et al.LESDEM study of 3-D motion of sand particles nearing bed in aeolian sand transport［J］.Acta Aerodynamica Sinica，2011，29（6）：784-788.

［3］ 李靜海，歐陽潔.顆粒流體復(fù)雜系統(tǒng)的多尺度模擬［M］.北京：科學出版社，2005：142-199.

［4］ CHIESA M，MATHIESEN V，MELHEIM J A，et al.Numerical simulation of particulate flow by the Eulerian-Lagrangian and the Eulerian-Eulerian approach with application to a fluidized bed［J］.Computers and Chemical Engineering，2005，29（2）：291-304.

［5］ 劉大有，董飛，賀大良.風沙二相流運動特點的分析［J］.地理學報，1996，51（5）：434-444.LIU Dayou，DONG Fei，HE Daliang.An analysis of the characters of blown sands［J］.Acta Geographica Sinica，1996，51（5）：434-444.

［6］ 顧兆林.風揚粉塵-近地層湍流與氣固兩相流 ［M］.北京：科學出版社，2010：79.

［7］ 顧正萌，郭烈錦.沙粒躍移運動的動理學模擬 ［J］.工程熱物理學報，2004，25（S1）：79-82.GU Zhengmeng，GUO Liejin.Simulation of sand saltation flow with kinetic theory ［J］.Journal of Engineering Thermophysics，2004，25（S1）：79-82.

［8］ 周芳，祁海鷹，由長福，等.有限空間風沙流動數(shù)值模擬及邊界條件問題 ［J］.清華大學學報：自然科學版，2004，44（8）：1079-1082.ZHOU Fang，QI Haiying，YOU Changfu，et al.Numerical simulation of wind-sand current and boundary conditions in a limited space［J］.Journal of Tsinghua University：Sci &Tech，2004，44（8）：1079-1082.

［9］ 武生智，任春勇.基于歐拉雙流體模型的風沙運動模擬 ［J］.蘭州大學學報：自然科學版，2012，48（1）：104-107.WU Shengzhi，REN Chunyong.Numerical simulation of wind blown sand based on the Eulerian model［J］.Journal of Lanzhou University： Natural Sciences，2012，48（1）：104-107.

［10］PIOMELLI U，BALARAS E.Wall-layer models for large eddy simulations ［J］.Annu Rev Fluid Mech，2002，34（1）：349-374.

［11］LU H L，GIDASPOW D.Hydrodynamics of binary fluidization in a riser：CFD simulation using two granular temperatures［J］.Chemical Engineering Science，2003，58（16）：3777-3792.

［12］FELICE R D.The voidage function for fluid-particle interaction systems［J］.International Journal of Multiphase Flow，1994，20（1）：153-159.

［13］曾卓雄.稠密兩相流動湍流模型及其應(yīng)用 ［M］.北京：機械工業(yè)出版社，2012：10-19.

［14］SYAMLAL M.MFIX documentation numerical technique，DOE／MC 31346-5824 ［R］.Washington DC，USA：DOE，1998：24.

［15］JOHNSON P C，JACKSON R.Frictional-collisional constitutive relations for granular materials，with application to plane shearing ［J］.Journal of Fluid Mechanics，1987，176（1）：67-93.

［16］ZHANG W，WANG Y，LEE S J.Simultaneous PIV and PTV measurements of wind and sand particle velocities［J］.Exp Fluids，2008，45（2）：241-256.

［17］LIU X P，DONG Z B.Experimental investigation of the concentration profile of a blowing sand cloud ［J］.Geomorphology，2004，60（3／4）：371-381.

［18］MATHIESEN V，SOLBERG T，HJERTAGER B H.An experimental and computational study of multiphase flow behavior in a circulating fluidized bed ［J］.International Journal of Multiphase Flow，2000，26（3）：387-419.

［19］KANG L Q，GUO L J.Wind tunnel experimental investigation of sand velocity in aeolian sand transport［J］.Geomorphology，2008，97（3／4）：438-450.

［20］HAFF P K，ANDERSON R S.Grain scale simulations of loose sedimentary beds：the example of grainbed impacts in aeolian saltation ［J］.Sedimentology，1993，40（2）：175-198.