近床面風(fēng)沙流的顆粒擬流體大渦模擬分析

周曉斯，王元，李志強(qiáng)

（西安交通大學(xué)流體機(jī)械及工程系，710049，西安）

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

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

連續(xù)性方程

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

氣相動(dòng)量方程

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

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

沙相動(dòng)量方程

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

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

針對(duì)顆粒擬溫度，建立如下方程

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

γs為顆粒脈動(dòng)能量引起的耗散率

Φs為氣固兩相之間的脈動(dòng)能量交換率

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

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

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

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

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

1.2 數(shù)值方法

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

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

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

2 計(jì)算模型及計(jì)算條件

2.1 計(jì)算模型

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

2.2 計(jì)算條件

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

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

計(jì)算域流向和展向均采用周期性邊界條件。在計(jì)算域法向上，沙床面處氣相速度采用無(wú)滑移條件，沙床面頂部采用對(duì)稱(chēng)邊界條件。計(jì)算時(shí)初始?jí)毫?chǎng)為0。借鑒文獻(xiàn)［2］，計(jì)算時(shí)間步長(zhǎng)為0.01ms。

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

顆粒擬溫度

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

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

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

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

3.1 來(lái)流風(fēng)場(chǎng)分析

圖3 凈風(fēng)場(chǎng)時(shí)均風(fēng)速廓線(xiàn)

3.2 沙相流場(chǎng)分析

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

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

2種曵力模型下躍移層沙相流向速度的統(tǒng)計(jì)平均值沿y的變化如圖5所示。為便于對(duì)比分析，圖中還給出了Liu等風(fēng)洞實(shí)驗(yàn)測(cè)得的沙相平均流向速度統(tǒng)計(jì)值的變化［17］。由圖5可見(jiàn)，躍移層內(nèi)，2種曵力模型在近床面的沙相速度均比較小，隨著y的增加，速度不斷提高。這是躍移沙粒不斷被來(lái)流風(fēng)場(chǎng)加速的緣故。在躍移中下層，沙相速度沿y變化的模擬值與實(shí)驗(yàn)值較為接近，在躍移上層二者偏差較大，模擬值大于實(shí)驗(yàn)值。文獻(xiàn)［17］的實(shí)驗(yàn)數(shù)據(jù)表明，相同自由來(lái)流風(fēng)速下，躍移沙相平均流向速度隨沙粒粒徑的增大而下降。在近床面，沙粒與床面以及沙粒之間的碰撞主導(dǎo)著沙粒運(yùn)動(dòng)；在離開(kāi)床面較高處，來(lái)流風(fēng)場(chǎng)成為沙粒運(yùn)動(dòng)的主導(dǎo)者，沙粒粒徑越小，沙粒與流體的跟隨性越強(qiáng)，沙粒易被加速，沙粒受風(fēng)加速的時(shí)間越長(zhǎng)，沙粒粒徑的影響導(dǎo)致的偏差凸顯。以上規(guī)律可以解釋本文模擬采用的ds=110μm單一沙粒粒徑下計(jì)算得到的在躍移層較高處的沙相平均流向速度，及Liu等在100μm＜ds≤200μm混合沙粒粒徑下得到的速度實(shí)驗(yàn)值出現(xiàn)較大偏差的原因。由圖5還可見(jiàn)，2種曵力模型下計(jì)算得到的速度沿y的變化達(dá)到高度吻合。這是因?yàn)橛?jì)算時(shí)風(fēng)沙流中的沙相稀疏，該狀態(tài)下2種曵力模型的系數(shù)相當(dāng)，所以2種模型的模擬結(jié)果達(dá)到一致。沙相平均流向速度沿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)計(jì)平均值沿y的變化如圖6所示。由圖6可見(jiàn)，與Liu等實(shí)驗(yàn)結(jié)果相比［17］，2種曵力模型的沙相濃度在躍移中下層的模擬值略高，在躍移上層偏低。這是因?yàn)楸疚哪M時(shí)在計(jì)算域流向和展向上均設(shè)定了周期性邊界條件，這2個(gè)方向上被夾帶流出區(qū)域的沙粒又由各個(gè)方向的入口返回，所以整個(gè)模擬過(guò)程中計(jì)算域平均沙相濃度保持不變。如果沙相濃度在躍移中下層被高估，則在躍移上層必然被低估。導(dǎo)致躍移上層出現(xiàn)偏差的另一種原因是，這2種曵力模型均低估了曵力系數(shù)，使得沙粒加速度偏小，沙粒弛豫時(shí)間加長(zhǎng)，大量沙粒不能被來(lái)流風(fēng)夾帶到高層，從而出現(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等的實(shí)驗(yàn)進(jìn)行對(duì)比可見(jiàn)，在相同來(lái)流風(fēng)速及合理粒徑條件下，擬流體模擬得到的沙相速度及濃度分布在宏觀統(tǒng)計(jì)值上均與實(shí)驗(yàn)值存在一定的偏差。為了實(shí)現(xiàn)氣沙兩相流場(chǎng)擬流體模型的準(zhǔn)確模擬，需要對(duì)風(fēng)沙流特性的擬流體氣固本構(gòu)封閉方程和相間作用模型進(jìn)行改進(jìn)。

定義沙相流向脈動(dòng)均方根速度（RMS）［18］

圖7 沙相uRMS沿y的變化

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

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

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

圖9 Felice曵力模型下e對(duì)uRMS的影響

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

4 結(jié) 論

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

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

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

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

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

［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］ 李志強(qiáng)，王元，王麗，等.風(fēng)沙流中近床面沙粒三維運(yùn)動(dòng)的LES-DEM分析 ［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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］ 李靜海，歐陽(yáng)潔.顆粒流體復(fù)雜系統(tǒng)的多尺度模擬［M］.北京：科學(xué)出版社，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］ 劉大有，董飛，賀大良.風(fēng)沙二相流運(yùn)動(dòng)特點(diǎn)的分析［J］.地理學(xué)報(bào)，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］ 顧兆林.風(fēng)揚(yáng)粉塵-近地層湍流與氣固兩相流 ［M］.北京：科學(xué)出版社，2010：79.

［7］ 顧正萌，郭烈錦.沙粒躍移運(yùn)動(dòng)的動(dòng)理學(xué)模擬 ［J］.工程熱物理學(xué)報(bào)，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］ 周芳，祁海鷹，由長(zhǎng)福，等.有限空間風(fēng)沙流動(dòng)數(shù)值模擬及邊界條件問(wèn)題 ［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，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］ 武生智，任春勇.基于歐拉雙流體模型的風(fēng)沙運(yùn)動(dòng)模擬 ［J］.蘭州大學(xué)學(xué)報(bào)：自然科學(xué)版，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］曾卓雄.稠密兩相流動(dòng)湍流模型及其應(yīng)用 ［M］.北京：機(jī)械工業(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.