風(fēng)沙流結(jié)構(gòu)的SPH模擬研究

楚花明，金阿芳，聞騰騰，梁嵐博

（新疆大學(xué)機(jī)械工程學(xué)院，新疆 烏魯木齊 830047）

1 引言

受人類(lèi)活動(dòng)的影響，自然原有環(huán)境遭到了極大破壞，土地沙漠化日益加劇，沙塵災(zāi)害成為世界面臨的嚴(yán)重自然災(zāi)害之一，極大的影響了人類(lèi)生存環(huán)境和生活質(zhì)量。在風(fēng)沙運(yùn)動(dòng)的研究中，輸沙量是風(fēng)沙流結(jié)構(gòu)的一個(gè)重要參數(shù)，它為制定有效的防沙治沙措施提供科學(xué)指導(dǎo)[1]，是各國(guó)學(xué)者研究的重點(diǎn)。目前對(duì)輸沙量沿高程上的分布規(guī)律不同學(xué)者有各自的研究結(jié)論還未形成共識(shí)。根據(jù)文獻(xiàn)[2]的研究，輸沙量與高度表現(xiàn)指數(shù)關(guān)系，即輸沙量隨高度的增加呈指數(shù)減小變化。他的結(jié)論被文獻(xiàn)[3]通過(guò)實(shí)驗(yàn)所證實(shí)并加以引用。然而在文獻(xiàn)[4]的研究中，他們發(fā)現(xiàn)輸沙量在高度上的分布并不嚴(yán)格遵守指數(shù)分布。當(dāng)前對(duì)風(fēng)沙流結(jié)構(gòu)的研究主要通過(guò)野外觀測(cè)、風(fēng)洞試驗(yàn)和數(shù)值模擬兩種手段進(jìn)行。目前國(guó)內(nèi)外關(guān)于輸沙量有三種表達(dá)式，國(guó)內(nèi)的主要結(jié)論還是指數(shù)函數(shù)。但受地表屬性的影響，近地層5cm內(nèi)可能發(fā)生明顯變化，但上層還是沒(méi)有脫離指數(shù)函數(shù)的變化規(guī)律。由于風(fēng)沙輸運(yùn)過(guò)程主要集中在近地層表面，沙粒濃度較高，通過(guò)風(fēng)洞實(shí)驗(yàn)觀察測(cè)量很難得到精準(zhǔn)的結(jié)果。文獻(xiàn)[5]的風(fēng)洞實(shí)驗(yàn)結(jié)果顯示，在近地層附近，輸沙量比指數(shù)分布數(shù)值要小。因此，輸沙量沿高度服從指數(shù)分布是在一定條件下得到的結(jié)果，并不是一個(gè)普適性結(jié)論。以往的研究結(jié)果表明，輸沙量在高程的分布結(jié)構(gòu)應(yīng)該與沙粒起跳時(shí)的速度分布存在關(guān)聯(lián)。因此，在研究風(fēng)沙流結(jié)構(gòu)時(shí)對(duì)沙粒起跳速度的分布規(guī)律的研究也必不可少。

現(xiàn)有對(duì)風(fēng)沙流運(yùn)動(dòng)的研究方法中主要包括：野外實(shí)地觀測(cè)，風(fēng)洞實(shí)驗(yàn)和數(shù)值模擬這三種方式[6]。隨著計(jì)算機(jī)硬件和軟件以及計(jì)算方法的發(fā)展，數(shù)值模擬方法因?yàn)榭梢员苊鈱?shí)地觀測(cè)的隨機(jī)性和不可復(fù)制性，以及相對(duì)于風(fēng)洞實(shí)驗(yàn)有成本較低和便于實(shí)現(xiàn)的優(yōu)點(diǎn)而在風(fēng)沙流運(yùn)動(dòng)的研究中越來(lái)越普及。光滑粒子流體動(dòng)力學(xué)（SPH）方法作為一種無(wú)網(wǎng)格拉格朗日粒子法，將待處理的問(wèn)題區(qū)域離散為一系列粒子，這些粒子既可以作為計(jì)算插值點(diǎn)，其身也攜帶有相關(guān)物理信息，特別適合處理水動(dòng)力學(xué)、爆炸沖擊以及天體運(yùn)動(dòng)等大變形問(wèn)題，應(yīng)用于風(fēng)沙流運(yùn)動(dòng)的研究上具有獨(dú)特的優(yōu)勢(shì)［7-8］。運(yùn)用SPH方法將氣體與沙粒均離散為SPH粒子，并建立氣固兩相之間的耦合關(guān)系，模擬了風(fēng)沙流運(yùn)動(dòng)過(guò)程并統(tǒng)計(jì)分析了沙粒的起跳速度分布，以及風(fēng)沙流達(dá)到動(dòng)態(tài)穩(wěn)定后的輸沙量沿高度的分布特點(diǎn)等風(fēng)沙流結(jié)構(gòu)的主要研究?jī)?nèi)容。

2 研究?jī)?nèi)容與方法

2.1 SPH方法介紹

SPH方法最早是由文獻(xiàn)［9］提出并應(yīng)用于模擬天體運(yùn)動(dòng)上的。它是將連續(xù)的物體離散為一系列相互聯(lián)系的SPH粒子，這些粒子攜帶有物體的諸如密度、硬度等物理信息，通過(guò)SPH控制方程來(lái)描述物體的宏觀性質(zhì)。其求解過(guò)程一般包括核近似和粒子近似兩個(gè)步驟，基本方程分別如下［10］。

2.1.1 核近似

對(duì)任意函數(shù)使用積分來(lái)近似表示，形如f(x)的函數(shù)，其積分表示式為：

式中：Ω—包含x的積分體積；x′—x影響域內(nèi)粒子；δ(x-x′)—Dirac函數(shù)，有以下性質(zhì)：

其中，使用核函數(shù)代替w(x-x′，h)，即狄拉克函數(shù)完成了核近似過(guò)程。如下式：

式中：w(x-x′，h)—核函數(shù)；h—光滑半徑。

2.1.2 粒子近似

通過(guò)對(duì)相鄰粒子的值進(jìn)行累加求和對(duì)函數(shù)的核近似表達(dá)式進(jìn)行近似，即產(chǎn)生目標(biāo)粒子處函數(shù)的粒子近似式。假設(shè)粒子的體積為ΔVj，則粒子的質(zhì)量為：mj=ΔVjρj（式中：ρj—粒子j的密度）。用ΔVj取代前述方程中粒子j處的無(wú)窮小體元dx′，則f(x)的連續(xù)SPH積分表示式為：

最終，粒子i處函數(shù)的粒子近似式可寫(xiě)為：

式中：WIJ=W(xi-xj，h)；N—粒子i的支持域內(nèi)的粒子總數(shù)。

2.2 風(fēng)沙運(yùn)動(dòng)模型

在風(fēng)沙運(yùn)動(dòng)中，氣流拖動(dòng)床面沙粒使其由靜止進(jìn)入運(yùn)動(dòng)狀態(tài)，而躍起的沙粒進(jìn)入空氣后獲得能量，在回落到床面時(shí)又擊濺起其他沙粒起跳，從而形成穩(wěn)定的運(yùn)動(dòng)狀態(tài)［11］。在這個(gè)過(guò)程中，風(fēng)沙流運(yùn)動(dòng)的一切能量都是由氣體提供的，而沙粒對(duì)氣體的運(yùn)動(dòng)具有反作用，改變氣體原先的運(yùn)動(dòng)狀態(tài)，從而形成了風(fēng)沙流這一復(fù)雜的、非線性且具有自組織性的兩相流系統(tǒng)［12］。

2.2.1 氣體控制方程

對(duì)于氣體，使用Navier-Stokes方程來(lái)描述其運(yùn)動(dòng)規(guī)律，其主要包括質(zhì)量守恒方程、動(dòng)量守恒方程公式。

式中：uβ—風(fēng)速在β方向上的分量；ρ—?dú)怏w密度。

式中：ν=μ/ρ—運(yùn)動(dòng)粘度；μ—?dú)怏w的動(dòng)力粘度；P和ρ—壓力和密度；u—?dú)怏w速度；α—方向分量。

式中：Rα—摩爾氣體常數(shù)，對(duì)于干燥氣體，Rα=287J/(kg·K)。

2.2.2 沙粒的運(yùn)動(dòng)方程

沙粒在運(yùn)動(dòng)過(guò)程中的受力情況，如圖1所示。從圖中可以看出沙粒在空中的受力情況較為復(fù)雜，所以將在沙粒運(yùn)動(dòng)過(guò)程中作用不明顯和目前尚無(wú)法量化的力，諸如薩夫曼力FS、Basset力FB以及電磁力忽略掉，則可以得到沙粒的運(yùn)動(dòng)方程為［13］：

圖1 沙粒受力分析Fig.1 Force Analysis of Sand

式中：FD—?dú)怏w拖曳力；FL—?dú)饬魃?；Fg—沙粒自身重力；FM—馬格納斯力。

2.2.3 沙粒的碰撞

沙粒之間的碰撞是風(fēng)沙運(yùn)動(dòng)過(guò)程中的一個(gè)重要機(jī)制，沙粒在落回床面以及在起跳進(jìn)入空中運(yùn)動(dòng)后有一定概率與其他沙粒發(fā)生碰撞。將發(fā)生碰撞的沙粒分別命名為i，j，k，i是一顆發(fā)生躍移后從空中降落的沙粒，如圖2所示。圖中：j—正在沙床上做蠕移運(yùn)動(dòng)的沙粒；k—沙床上的一顆隨機(jī)沙粒。假設(shè)，沙粒i以水平夾角為α的速度Vi與以速度Vj運(yùn)動(dòng)的沙粒j發(fā)生碰撞，而沙粒j又與沙床上靜止的沙粒k發(fā)生碰撞。碰撞點(diǎn)如圖中標(biāo)識(shí)為1，2。沙粒i與沙粒j碰撞時(shí)兩沙粒的質(zhì)心連線和水平方向的夾角記為β，沙粒j與沙粒k的質(zhì)心連線與豎直方向的夾角記為γ。

圖2 沙粒與沙床以及與空中沙粒的碰撞Fig.2 The Collision of Sand-Bed and Sand-Sand

不同于其他模擬方法單獨(dú)設(shè)定碰撞概率來(lái)模擬碰撞，SPH方法本身就會(huì)在每一時(shí)間步自動(dòng)搜索目標(biāo)粒子支持域范圍內(nèi)的沙粒，對(duì)支持域內(nèi)判定為沙粒的粒子會(huì)根部光滑長(zhǎng)度和兩粒子之間的距離來(lái)判斷是否發(fā)生碰撞，對(duì)符合碰撞條件的使用動(dòng)量定理來(lái)計(jì)算出碰撞后兩粒子的速度。

2.3 計(jì)算模型及參數(shù)設(shè)置

主要是基于沙粒與氣體粒子之間的受力機(jī)制來(lái)模擬沙粒的運(yùn)動(dòng)過(guò)程并分析風(fēng)沙流結(jié)構(gòu)特點(diǎn)，分別模擬了在摩阻風(fēng)速為0.1（m/s）和0.5（m/s）時(shí)，0.1mm、0.25mm兩種粒徑沙粒的運(yùn)動(dòng)過(guò)程。模擬過(guò)程中涉及的各類(lèi)參數(shù)，如表1所示。

表1 算例中氣體與沙粒的參數(shù)列表Tab.1 Parameter List of Gas and S and in Numerical Example

在模擬區(qū)域的上下邊界各設(shè)置一層虛粒子防止粒子穿透，模擬區(qū)域的左右邊界采取周期性邊界條件模擬粒子的運(yùn)動(dòng)過(guò)程，如圖3所示。初始來(lái)流風(fēng)速遵循式（11）所示的對(duì)數(shù)分布[14]，如圖4所示。

圖3 計(jì)算區(qū)域初始狀態(tài)示意圖Fig.3 Initial State of Calculated Area

圖4 風(fēng)速廓線圖Fig.4 Wind Speed Profile

式中：uy—高度為y處的風(fēng)速；k—馮·卡曼系數(shù)；Y0—沙床面粗糙度。

3 結(jié)果與分析

3.1 模型驗(yàn)證

我們選取了6顆具有典型躍移沙粒運(yùn)動(dòng)特征的沙粒，所選粒子位置，如圖5所示。把它們?cè)诟鱾€(gè)時(shí)間步的位置依次連接得到它們的運(yùn)動(dòng)軌跡，其中1、4、6號(hào)粒子為直徑0.1mm的沙粒，2、3、5號(hào)粒子為直徑0.25mm的沙粒，如圖6所示。從圖中可以看出，沙粒的躍移軌跡具有非常明顯的拋物線特性，這驗(yàn)證了我們建立的沙粒受力模型是正確的。沙粒的躍移高度和長(zhǎng)度會(huì)因沙粒的起跳速度和角度不同而不同：圖中1號(hào)和5號(hào)沙粒在躍移結(jié)束后落回床面不再運(yùn)動(dòng)；2號(hào)沙粒在落回床面后因與床面上的沙粒發(fā)生碰撞而反彈，反彈后會(huì)做非常短距離的躍移；3號(hào)、4號(hào)和6號(hào)沙粒在落回床面后以較大的反彈速度和角度重新彈回氣流中，其中6號(hào)沙粒的反彈速度很大；可見(jiàn)，大多數(shù)的躍移沙粒在落回床面時(shí)都會(huì)因與其他沙粒的碰撞而反彈，但反彈的速度和角度根據(jù)每個(gè)沙粒的初始速度和位置不同有較大差異。

圖5 所選粒子位置示意圖Fig.5 Position of Selected Particle

圖6 沙粒的躍移軌跡Fig.6 Jump Trajectory of S and Particles

3.2 沙粒的起跳速度分布

在風(fēng)沙流結(jié)構(gòu)的研究中，沙粒的起跳速度分布將影響到輸沙量沿高度的分布。我們模擬了沙粒在沙床表面被風(fēng)吹起，在空中運(yùn)動(dòng)后落回床面并擊起其它沙粒起跳的整個(gè)過(guò)程。隨機(jī)選取了100顆起跳沙粒，分別對(duì)它們起跳時(shí)的水平速度和垂直速度進(jìn)行統(tǒng)計(jì)分析，如圖7、圖8所示。從圖7（a）和圖7（b）可以看出，沙粒的水平起跳速度分布具有明顯的單峰特性，大部分的水平起跳速度分布在（0.1～0.4）m/s之間；同等粒徑的沙粒在來(lái)流風(fēng)速較大時(shí)，其水平起跳速度較大；0.1mm和0.25mm兩種粒徑沙粒的起跳速度分布區(qū)別不大。圖8（a）和圖8（b）顯示了躍移沙粒垂直起跳速度的概率密度分布?？梢?jiàn)，兩種粒徑的沙粒其躍移的垂直起跳速度呈指數(shù)分布，隨著垂直速度的增大起跳數(shù)量呈指數(shù)下降。相對(duì)而言，小粒徑沙粒的垂直起跳速度要大于同等摩阻風(fēng)速條件下的大粒徑沙粒，同一粒徑的沙粒在兩種摩阻風(fēng)速下的垂直起跳速度分布基本一致，說(shuō)明沙粒粒徑對(duì)垂直起跳的速度影響更大。

圖7 沙粒起跳水平速度概率分布Fig.7 Distribution of Horizontal of Take-Off Sand

圖8 沙粒起跳垂直速度概率分布Fig.8 Distribution of Vertical Velocity of Take-Off Sand

不同粒徑沙粒在不同摩阻風(fēng)速下的起跳速度概率分布，如圖9所示?？梢钥闯觯S著起跳速度的增加，沙粒的起跳數(shù)量經(jīng)歷了一個(gè)先增加達(dá)到峰值后減小的過(guò)程，單峰向左偏移，這符合gamma分布。這與Anderson和Hallet等人的研究相吻合[13]。同時(shí)，在同一粒徑條件下，摩阻風(fēng)速為0.5（ms-1）時(shí)，其起跳速度分布函數(shù)概率密度相對(duì)于摩阻風(fēng)速為0.1m/s時(shí)要向右分散。

圖9 躍移沙粒起跳速度的分布Fig.9 Take-off Velocity Distribution of Leaping S and Particle

3.3 單寬輸沙率

當(dāng)沙床面上起跳的沙粒數(shù)與落回床面的沙粒數(shù)相差不大時(shí)，風(fēng)沙運(yùn)動(dòng)進(jìn)入了穩(wěn)態(tài)。兩種粒徑沙粒在不同摩阻風(fēng)速下，單寬輸沙率隨時(shí)間的變化關(guān)系，如圖10所示。

圖10 單寬輸沙率隨時(shí)間步發(fā)展的過(guò)程Fig.10 The Process of Single Width S and Transport with Time

從圖中可知，在相同粒徑下，摩阻風(fēng)速為0.5（ms-1）的單寬輸沙率要顯著大于摩阻風(fēng)速為0.1（ms-1）時(shí)，而在相同的摩阻風(fēng)速下，粒徑較小沙粒的單寬輸沙率要大于粒徑較大的沙粒；不同粒徑沙粒的單寬輸沙率都隨著時(shí)間步的增加而先增加后稍有減小并趨于穩(wěn)定，在一定時(shí)間步后，單寬輸沙率穩(wěn)定在一數(shù)值上下。當(dāng)風(fēng)沙運(yùn)動(dòng)處于動(dòng)態(tài)平衡時(shí)，0.1mm粒徑沙粒在摩阻風(fēng)速0.1（ms-1）和0.5（ms-1）下的單寬輸沙率分別為：0.015（kgm-1s-1）、0.0461（kgm-1s-1），0.25mm粒徑沙粒分別為0.0074（kgm-1s-1），0.038（kgm-1s-1）。由圖中可以看到，0.1mm粒徑沙粒的單寬輸沙率在兩種摩阻風(fēng)速下的增長(zhǎng)速度都高于0.25mm粒徑沙粒，在0.5（ms-1）摩阻風(fēng)速下更快達(dá)到動(dòng)態(tài)平衡。但在0.1（ms-1）摩阻風(fēng)速下，0.25mm粒徑沙粒要先達(dá)到動(dòng)態(tài)平衡，這與風(fēng)速較小且0.25mm粒徑沙粒起跳并不充分有關(guān)，其在起跳沙粒數(shù)量較少的情況下達(dá)到動(dòng)態(tài)平衡。

3.4 輸沙量

當(dāng)沙床上沙粒的吹蝕和堆積數(shù)量相當(dāng)時(shí)，風(fēng)沙運(yùn)動(dòng)處于動(dòng)態(tài)平衡，此時(shí)輸沙量的大小可由該高度處沙粒的數(shù)量來(lái)表示，因此，輸沙量沿高度的變化即為在單位時(shí)間內(nèi)單位面積上沙粒濃度沿高度的變化。由上文得出沙粒的起跳速度服從gamma分布，可見(jiàn)，沙粒以不同的速度進(jìn)入到空中，其所能達(dá)到的最大高度也不同。文中統(tǒng)計(jì)了在風(fēng)沙流穩(wěn)定后，輸沙量沿高度的變化。隨著高度的增加，輸沙量先增加而后減小呈指數(shù)衰減，在接近沙床面的上方附近存在一個(gè)最大值，如圖10所示。

這與文獻(xiàn)[14]所做的模擬結(jié)果趨勢(shì)上相同。由圖11（a）和圖11（b）可以看出，在沙粒粒徑為0.1mm和0.25mm兩種情況下，摩阻風(fēng)速為0.5（ms-1）時(shí)，其輸沙量達(dá)到最大值的高度要高于摩阻風(fēng)速為0.1（ms-1）時(shí)，輸沙量也要遠(yuǎn)大于摩阻風(fēng)速為0.1（ms-1）時(shí)，這是因?yàn)殡S著風(fēng)速的增加沙粒起跳時(shí)攜帶的能量也增加，所以輸沙量的最大值會(huì)隨著風(fēng)速的增大而增大，即“象鼻效應(yīng)”。

圖11 輸沙量沿高度的分布Fig.11 Distribution of S and Transport Along Height

4 討論

（1）影響沙粒起跳速度的因素有很多，文中只考慮了0.1mm和0.25mm兩種粒徑沙粒在0.1（ms-1）和0.5（ms-1）兩種摩阻風(fēng)速下的運(yùn)動(dòng)過(guò)程。由起跳速度的分布情況可知，在文中條件下沙粒的粒徑對(duì)起跳的速度影響更大，這是因?yàn)樵谏沉Ａ綖?.1mm時(shí)，給出的兩種摩阻風(fēng)速下都能夠充分起跳；沙粒粒徑為0.25mm時(shí)，兩種摩阻風(fēng)速下的起跳都不夠充分。

（2）輸沙量沿高度的分布和沙粒的起跳速度分布具有很大的相似性，沙粒的起跳速度具有一單峰值，并偏移在低速度區(qū)域內(nèi)，輸沙量則沿著高度先增加而后指數(shù)減少，也存在一個(gè)輸沙量峰值，同樣偏移在距離地面較近的高度，該高度值與沙粒的起跳速度分布有關(guān)，在沙粒的起跳速度分布中存在一個(gè)速度區(qū)間，在該速度區(qū)間內(nèi)，起跳的沙粒數(shù)量最多，以該速度起跳的沙粒所能達(dá)到的高度即為輸沙量達(dá)到最大值時(shí)的高度。

（3）不同于其他參數(shù)化建模的模擬方式，文中基于SPH方法建立的風(fēng)沙流模型是從沙粒及其氣體粒子的受力情況入手，將沙粒的受力情況量化，選取沙粒在運(yùn)動(dòng)過(guò)程中一直起主要作用的力，忽略掉一些次要的作用力，建立風(fēng)沙兩項(xiàng)間的耦合機(jī)制以及沙粒之間的碰撞機(jī)制，自然得到沙粒的起跳速度及其它一系列運(yùn)動(dòng)情況，更容易探究沙粒的受力情況，符合物理本質(zhì)。

5 結(jié)論

（1）不同粒徑的沙粒，其風(fēng)沙運(yùn)動(dòng)演變過(guò)程基本一致，其輸沙量都隨著摩阻風(fēng)速的增加而增加，隨著時(shí)間的增加而先呈指數(shù)形式的增長(zhǎng)而后略有減小，最后趨于穩(wěn)定。粒徑大的沙粒其輸沙量達(dá)到穩(wěn)定狀態(tài)需要的時(shí)間略長(zhǎng)。

（2）通過(guò)SPH方法對(duì)沙粒運(yùn)動(dòng)過(guò)程進(jìn)行模擬，發(fā)現(xiàn)其起跳速度呈gamma分布，輸沙量沿高度先增加而后呈指數(shù)減小。摩阻風(fēng)速較大時(shí)，輸沙量的達(dá)到的峰值以及該峰值對(duì)應(yīng)的高度都大于摩阻風(fēng)速較低時(shí)的相應(yīng)值。

（3）輸沙量沿高度的分布和沙粒起跳的速度分布有內(nèi)在關(guān)聯(lián)，輸沙量其飽和層的高度和沙粒起跳速度分布的峰值所對(duì)應(yīng)的速度有關(guān)，以該速度起跳所達(dá)到的高度即為飽和層高度。

（4）SPH方法基于自身特點(diǎn)非常適合風(fēng)沙兩相流的模擬。課題組首次運(yùn)用SPH方法，建立了風(fēng)沙兩相流的運(yùn)動(dòng)模型，模擬結(jié)果與理論吻合較好，在風(fēng)沙運(yùn)動(dòng)的機(jī)理研究中具有很大應(yīng)用前景。