風沙流結構的SPH模擬研究

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

（新疆大學機械工程學院，新疆 烏魯木齊 830047）

1 引言

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

現有對風沙流運動的研究方法中主要包括：野外實地觀測，風洞實驗和數值模擬這三種方式[6]。隨著計算機硬件和軟件以及計算方法的發展，數值模擬方法因為可以避免實地觀測的隨機性和不可復制性，以及相對于風洞實驗有成本較低和便于實現的優點而在風沙流運動的研究中越來越普及。光滑粒子流體動力學（SPH）方法作為一種無網格拉格朗日粒子法，將待處理的問題區域離散為一系列粒子，這些粒子既可以作為計算插值點，其身也攜帶有相關物理信息，特別適合處理水動力學、爆炸沖擊以及天體運動等大變形問題，應用于風沙流運動的研究上具有獨特的優勢［7-8］。運用SPH方法將氣體與沙粒均離散為SPH粒子，并建立氣固兩相之間的耦合關系，模擬了風沙流運動過程并統計分析了沙粒的起跳速度分布，以及風沙流達到動態穩定后的輸沙量沿高度的分布特點等風沙流結構的主要研究內容。

2 研究內容與方法

2.1 SPH方法介紹

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

2.1.1 核近似

對任意函數使用積分來近似表示，形如f(x)的函數，其積分表示式為：

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

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

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

2.1.2 粒子近似

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

最終，粒子i處函數的粒子近似式可寫為：

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

2.2 風沙運動模型

在風沙運動中，氣流拖動床面沙粒使其由靜止進入運動狀態，而躍起的沙粒進入空氣后獲得能量，在回落到床面時又擊濺起其他沙粒起跳，從而形成穩定的運動狀態［11］。在這個過程中，風沙流運動的一切能量都是由氣體提供的，而沙粒對氣體的運動具有反作用，改變氣體原先的運動狀態，從而形成了風沙流這一復雜的、非線性且具有自組織性的兩相流系統［12］。

2.2.1 氣體控制方程

對于氣體，使用Navier-Stokes方程來描述其運動規律，其主要包括質量守恒方程、動量守恒方程公式。

式中：uβ—風速在β方向上的分量；ρ—氣體密度。

式中：ν=μ/ρ—運動粘度；μ—氣體的動力粘度；P和ρ—壓力和密度；u—氣體速度；α—方向分量。

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

2.2.2 沙粒的運動方程

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

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

式中：FD—氣體拖曳力；FL—氣流升力；Fg—沙粒自身重力；FM—馬格納斯力。

2.2.3 沙粒的碰撞

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

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

不同于其他模擬方法單獨設定碰撞概率來模擬碰撞，SPH方法本身就會在每一時間步自動搜索目標粒子支持域范圍內的沙粒，對支持域內判定為沙粒的粒子會根部光滑長度和兩粒子之間的距離來判斷是否發生碰撞，對符合碰撞條件的使用動量定理來計算出碰撞后兩粒子的速度。

2.3 計算模型及參數設置

主要是基于沙粒與氣體粒子之間的受力機制來模擬沙粒的運動過程并分析風沙流結構特點，分別模擬了在摩阻風速為0.1（m/s）和0.5（m/s）時，0.1mm、0.25mm兩種粒徑沙粒的運動過程。模擬過程中涉及的各類參數，如表1所示。

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

在模擬區域的上下邊界各設置一層虛粒子防止粒子穿透，模擬區域的左右邊界采取周期性邊界條件模擬粒子的運動過程，如圖3所示。初始來流風速遵循式（11）所示的對數分布[14]，如圖4所示。

圖3 計算區域初始狀態示意圖Fig.3 Initial State of Calculated Area

圖4 風速廓線圖Fig.4 Wind Speed Profile

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

3 結果與分析

3.1 模型驗證

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

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

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

3.2 沙粒的起跳速度分布

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

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

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

不同粒徑沙粒在不同摩阻風速下的起跳速度概率分布，如圖9所示。可以看出，隨著起跳速度的增加，沙粒的起跳數量經歷了一個先增加達到峰值后減小的過程，單峰向左偏移，這符合gamma分布。這與Anderson和Hallet等人的研究相吻合[13]。同時，在同一粒徑條件下，摩阻風速為0.5（ms-1）時，其起跳速度分布函數概率密度相對于摩阻風速為0.1m/s時要向右分散。

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

3.3 單寬輸沙率

當沙床面上起跳的沙粒數與落回床面的沙粒數相差不大時，風沙運動進入了穩態。兩種粒徑沙粒在不同摩阻風速下，單寬輸沙率隨時間的變化關系，如圖10所示。

圖10 單寬輸沙率隨時間步發展的過程Fig.10 The Process of Single Width S and Transport with Time

從圖中可知，在相同粒徑下，摩阻風速為0.5（ms-1）的單寬輸沙率要顯著大于摩阻風速為0.1（ms-1）時，而在相同的摩阻風速下，粒徑較小沙粒的單寬輸沙率要大于粒徑較大的沙粒；不同粒徑沙粒的單寬輸沙率都隨著時間步的增加而先增加后稍有減小并趨于穩定，在一定時間步后，單寬輸沙率穩定在一數值上下。當風沙運動處于動態平衡時，0.1mm粒徑沙粒在摩阻風速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粒徑沙粒的單寬輸沙率在兩種摩阻風速下的增長速度都高于0.25mm粒徑沙粒，在0.5（ms-1）摩阻風速下更快達到動態平衡。但在0.1（ms-1）摩阻風速下，0.25mm粒徑沙粒要先達到動態平衡，這與風速較小且0.25mm粒徑沙粒起跳并不充分有關，其在起跳沙粒數量較少的情況下達到動態平衡。

3.4 輸沙量

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

這與文獻[14]所做的模擬結果趨勢上相同。由圖11（a）和圖11（b）可以看出，在沙粒粒徑為0.1mm和0.25mm兩種情況下，摩阻風速為0.5（ms-1）時，其輸沙量達到最大值的高度要高于摩阻風速為0.1（ms-1）時，輸沙量也要遠大于摩阻風速為0.1（ms-1）時，這是因為隨著風速的增加沙粒起跳時攜帶的能量也增加，所以輸沙量的最大值會隨著風速的增大而增大，即“象鼻效應”。

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

4 討論

（1）影響沙粒起跳速度的因素有很多，文中只考慮了0.1mm和0.25mm兩種粒徑沙粒在0.1（ms-1）和0.5（ms-1）兩種摩阻風速下的運動過程。由起跳速度的分布情況可知，在文中條件下沙粒的粒徑對起跳的速度影響更大，這是因為在沙粒粒徑為0.1mm時，給出的兩種摩阻風速下都能夠充分起跳；沙粒粒徑為0.25mm時，兩種摩阻風速下的起跳都不夠充分。

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

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

5 結論

（1）不同粒徑的沙粒，其風沙運動演變過程基本一致，其輸沙量都隨著摩阻風速的增加而增加，隨著時間的增加而先呈指數形式的增長而后略有減小，最后趨于穩定。粒徑大的沙粒其輸沙量達到穩定狀態需要的時間略長。

（2）通過SPH方法對沙粒運動過程進行模擬，發現其起跳速度呈gamma分布，輸沙量沿高度先增加而后呈指數減小。摩阻風速較大時，輸沙量的達到的峰值以及該峰值對應的高度都大于摩阻風速較低時的相應值。

（3）輸沙量沿高度的分布和沙粒起跳的速度分布有內在關聯，輸沙量其飽和層的高度和沙粒起跳速度分布的峰值所對應的速度有關，以該速度起跳所達到的高度即為飽和層高度。

（4）SPH方法基于自身特點非常適合風沙兩相流的模擬。課題組首次運用SPH方法，建立了風沙兩相流的運動模型，模擬結果與理論吻合較好，在風沙運動的機理研究中具有很大應用前景。