SPH方法在風沙流研究中的應用進展

趙杰，金阿芳，楚花明

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

光滑粒子流體動力學(Smoothed Particle Hydrodynamics)又被稱為SPH方法，是在1977年最早由Gingold和Monaghan等人提出來并加以運用的[1,2]。該方法最初是應用在天體物理中以解決流體質團在無窮大三維空間中任意流動的問題。隨著其他研究人員和學者的應用以及改進，SPH方法快速發展，并開始推廣應用于流體力學方面的研究。SPH作為一種數值模擬方法同時具有無網格性質和粒子性質，特別適合于處理大變形問題、多相流問題以及追蹤單個粒子運動[3]，在潰壩、水下爆炸、飛沫方面的問題研究中都取得了成果。

另一方面，為了研究風沙流的運動機理，防治風沙災害，對風沙運動的研究在20世紀30年代初便引起了當時西方學者的關注。英國學者拜格諾(Bagnold)在1941發表了自己第一篇關于風沙地貌學的著作，標志著風沙物理學的開端[4]。此后，關于風沙運動的研究進一步發展，逐漸發展到現有的3種研究方法：野外觀測、風洞實驗和數值模擬。而對風沙流的研究也著眼于宏觀和微觀兩個方面：宏觀上注重對風沙流整體結構的研究，微觀上則對單顆沙粒的運動規律進行研究，其中，躍移沙粒占整個輸移沙粒總量的75%，成為連接宏觀與微觀研究之間的橋梁[5]。在現有的研究中，尚未對風沙運動的內在機理形成統一認識，對風沙流的定量描述也多是基于觀測數據的半經驗公式。所以，使用SPH方法對風沙運動進行數值模擬方法將有助于研究風沙流的運動規律并幫助完善其基本理論，同時也拓寬了SPH方法的應用領域，促進方法本身的成熟。

1 SPH模擬風沙流

風沙流是一種典型的氣固兩相流，其中包括固體和氣體兩種相互作用的物質。SPH方法作為一種流體擬顆粒模型，在將沙粒相作為離散相來處理同時，將氣體也離散為“顆粒”性質的流體微團，通過兩種不同物質顆粒之間的相互作用，來實現風沙流物理特性的模擬[6]。由于在模擬過程中不需要建立網格，同時又能攜帶單個粒子的信息，所以特別適合處理風沙流這種大變形的物理運動。

金阿芳等人從氣固兩相流的力學觀點出發，通過沙粒相和氣體相之間的碰撞建立了兩相之間的耦合關系，建立了基于光滑粒子流體動力學方法的風沙運動控制方程[5]。在他們的工作中，編程實現了SPH算法對風沙流運動的模擬，提出了使用SPH方法對風沙流模擬時的一些關鍵參數的選取規則，在處理邊界條件時，他們使用虛粒子來阻止粒子的穿透，為今后在對風沙流動進行數值模擬時，確定邊界條件給出了一定的參考。同時，他們通過仿真模擬了躍移沙粒運動狀態，確定了其運動參數和幾何參數，其中運動參數包括起跳速度、沉降速度、碰撞速度等。幾何參數包括躍移長度和高度、起跳角和碰撞角等，并對計算結果進行了分析。 牟陽陽等人使用SPH方法建立了傾斜沙床面下風沙運動的物理模型，他們模擬分析了在不同時間段沙粒在空間中的分布特點；統計了在沙床斜面上不同位置沙粒的起跳高度和距離，得出沙粒的平均躍移距離是單個沙波紋長度的1.16倍[7]。在逯博等人的研究中，對氣體相使用N-S方程，對沙粒相使用牛頓定律進行處理，建立了兩相的控制方程，他們使用SPH方法對沙粒的躍移進行了模擬，得到了沙粒的典型躍移軌跡，并在模擬中加入了風的加載和卸載模塊，統計了在風的加載和卸載情況下沙粒的軌跡和法向速度的變化[8]。為了對沙粒的躍移參數進行總體性地描述，金阿芳等人在模擬中隨機選取了1 000個沙粒進行分析，得到了它們的起跳角和起跳速度的概率密度分布。他們發現沙粒的躍移起跳角度的概率密度函數符合指數分布，起跳速度的概率密度函數符合對數正態分布；并在模擬中，他們對沙粒相和氣體相分別選取了不同的光滑長度，得到了更好的耦合效果[9]。

陳福振等人則將SPH方法和傳統的有限體積法(FVM)結合起來，建立了SPH-FVM耦合方法來模擬風沙運動[10]。在他們的研究中，對連續氣體相采用FVM，將沙粒相視為擬流體采用SPH方法處理，氣體相和沙粒相之間通過拖曳力、壓力梯度和濃度等參數進行耦合。同時，他們對SPH方法進行了改造，SPH粒子代表顆粒群，包含沙粒群諸如密度、質量、體積、粒徑分布和體積分數等物理信息，由此減小了問題規模，提高了計算效率。陳福振等人在SPH邊界條件的設置上使用了罰函數方法施加邊界力，對氣體相使用無滑移邊界條件。在對風沙流的模擬中，他們得到了沙粒的運動軌跡、沙粒在氣體中的水平速度分布以及氣體在沙粒的反作用下速度的變化，驗證了方法的有效性。

2 模型和參數選取

在已有的基于SPH方法建立的風沙運動的物理模型中，因為編程實現難度以及研究內容的側重點不同，都對模型做了一定程度的簡化。同時，因為現在尚未有對風沙運動的內在機理形成統一的共識，在處理物理模型時，使用的方法也不同。例如，在現有基于SPH方法的風沙流研究中都未考慮沙粒的帶電模型。在構建沙床面時，雖然大多數沙床面由混合粒徑的沙粒構成，但不同粒徑的沙粒在沙床面的隨機分布并不完全符合不同其在實際沙床上的分布規律。在研究沙粒的運動過程中，大多數模型都忽略了沙粒自身的旋轉和沙粒的蠕移，這將影響到沙粒的躍移高度以及當沙粒落回床面時，其擊濺起沙粒的起跳角度，同時，對沙床面的形態也會造成影響。對來流風的處理上，大多簡化為均勻流或對數流，對紊流的影響并未加以考慮。

當前，限于建立的物理模型規模，研究的主要內容仍是單顆沙粒運動的微觀特征，對大尺度的風沙運動特性如風沙流的輸移特征、沙丘的形成過程等還有待進一步的研究。在使用SPH方法建模時，對方法中一些重要參數的選取，如：核函數、光滑長度、時間步長等大都根據已有的經驗，或考慮到計算效率予以選擇，這些都將影響到模擬結果的精度[11]。

3 邊界條件

在SPH中，對邊界條件的處理仍然沒有一致的結論。在風沙運動的模擬中，金阿芳等人在他們的工作中通常在沙床面下方設定一個固定邊界，防止粒子向下穿透，氣體上方設定自由邊界或固定邊界，計算域的左右兩端往往是周期邊界，即運動出計算域的粒子將攜帶已有的物理信息從另一側重新進入[5]。自從Monaghan于1994年首次在自由表面流動中使用了排斥邊界后[12]，后來的學者又提出了一系列方法，這些方法可以歸為3類：虛粒子法、排斥函數和邊界積分。其中，因為虛粒子法和排斥函數的應用條件相對較為簡單，其使用較為普遍。虛粒子法是在邊界外排布一組虛粒子，參與問題域內粒子的粒子近似計算中，使得原問題域內粒子的支持域不被截斷，提高了近似精度。這種方法雖然在處理邊界條件時十分有效，但同時它無法完全阻止粒子穿透邊界，或會產生非物理分離。鏡像粒子是另外一種虛粒子方法，它通過將粒子沿著邊界生成的鏡像來防止粒子穿透，但這種方法只適用于簡單幾何形狀的邊界。排斥函數的方法則用另外一種方式使用虛粒子，它在邊界處布置一排虛粒子，對接近虛粒子的粒子產生排斥力，排斥力大小和兩種粒子之間的距離有關。Monaghan和Kajtar隨后對這種方法進行了改進，解決了當粒子平行于邊界運動時受力不穩定的問題[13]。但是這種方法的缺陷在于，靠近于邊界的粒子，其支持域被截斷，無法得到很精確的結果。

為了修正上述在邊界處存在的缺陷、提高邊界精度，不斷有學者提出對SPH方法的改進形式。如Dilts等人基于伽遼金近似構造的移動最小二乘光滑粒子法(MLSPH)[14]， Chen建立的正則化形式的修正光滑粒子法 ( CSPH)[15]， Liu使用泰勒展開在非連續區域兩端分段構造的非連續光滑粒子法 ( DSPH)[16]，Liu對核函數采用校正函數進行修正的再生核粒子方法 ( RKPM)[16]。總體而言，當前在發展SPH新的邊界處理方法方面取得了一定進展，改善了邊界出的計算精度。但是在邊界處理上仍存有許多未被發現的問題，有待進一步解決。

4 計算效率

在現有的基于SPH方法對風沙流的模擬中，大多數為二維模型。雖然因為SPH程序的特點，由二維環境擴展到三維并不困難，但粒子數的增加在提高結果精度的同時也將大幅增加計算量。為兼顧精度和計算效率，目前通常忽略風沙流在平行于地面、垂直于氣流方向上的變化，將風沙運動模型簡化為二維情況。在SPH程序的計算用時中，絕大部分時間是用于其粒子搜索過程。目前在風沙運動模擬中較常使用的全配對粒子搜索法最為簡單，其基本原理是：在SPH方法中若要得到粒子i的物理信息，首先要得到其支持域內所有粒子的信息，對其疊加求和后加權平均得到粒子i的近似解。所以在每個時間步的開始都需要在粒子i的支持域內進行粒子搜索，在更新了粒子i的信息后，再進行下一個時間步。當粒子的個數為n時，則每個粒子需要搜索n-1次，則每一步最終的計算量為 O( n2)。因此，在風沙流這種粒子數眾多的問題中，計算速度較慢。此外，其他的粒子搜索法還有鏈表搜索法，但當問題域在運算過程中發生大的變形時，這種方法仍然效率較低[17]。

近來，基于圖形處理單元(GPU)的開源代碼快速發展，逐漸應用到數值計算當中。因為SPH非常適合使用GPU進行數據并行處理，已經發展出來大量基于GPU的開源代碼，例如Hérault等人提出的GPUSPH[18]，Crespo等人的DualSPHysics[19]。它們可以在個人計算機上只用數小時便模擬上百萬個粒子，這對之前所需的幾個月時間是一個巨大的飛躍，也為SPH方法應用于工程實踐提供了可能。

5 總結

相較于傳統的網格方法，SPH以其無網格、具有粒子性質的特點使得其在應用于諸如風沙流等大變形問題上具有優勢。目前，SPH方法在計算穩定性和可行性上都取得了一定的進展。但若要將基于SPH的風沙流模擬應用到工程實踐中，仍需要解決以下問題：

5.1 對風沙運動的模擬還有賴于風沙物理學的進一步發展。當前風沙物理學仍處于發展過程，尚未形成完整成熟的理論。在使用控制方程對風沙兩相流進行數學描述時，還無法將影響因素考慮完全，如還未考慮能量耗散等問題。由此，使用SPH方法對風沙運動進行模擬時還需要風沙物理學理論的發展和完善。

5.2 使用SPH方法計算效率還需要進一步提高。由于SPH方法的計算量十分巨大，用SPH方法模擬風沙流還僅限于二維情況，對三維風沙運動的模擬還未出現。當前已有國內外學者針對GPU并行算法展開了研究并取得了巨大成果，未來將基于GPU的并行運算引入到風沙運動的模擬中，將進一步擴大模擬規模，提高模擬精度。

5.3 SPH方法還需要進一步完善。當前，SPH方法在某些領域的模擬中取得了較好的結果，但其還應該繼續完善以下4個方面：數值穩定性、收斂性、邊界條件以及自適應性。隨著SPH方法不斷地發展和完善，使用SPH方法對風沙運動進行模擬的精度將更加可靠。