趙杰，金阿芳，楚花明

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

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

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

1 SPH模擬風(fēng)沙流

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

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

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

2 模型和參數(shù)選取

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

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

3 邊界條件

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

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

4 計算效率

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

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

5 總結(jié)

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

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

5.2 使用SPH方法計算效率還需要進(jìn)一步提高。由于SPH方法的計算量十分巨大，用SPH方法模擬風(fēng)沙流還僅限于二維情況，對三維風(fēng)沙運動的模擬還未出現(xiàn)。當(dāng)前已有國內(nèi)外學(xué)者針對GPU并行算法展開了研究并取得了巨大成果，未來將基于GPU的并行運算引入到風(fēng)沙運動的模擬中，將進(jìn)一步擴(kuò)大模擬規(guī)模，提高模擬精度。

5.3 SPH方法還需要進(jìn)一步完善。當(dāng)前，SPH方法在某些領(lǐng)域的模擬中取得了較好的結(jié)果，但其還應(yīng)該繼續(xù)完善以下4個方面：數(shù)值穩(wěn)定性、收斂性、邊界條件以及自適應(yīng)性。隨著SPH方法不斷地發(fā)展和完善，使用SPH方法對風(fēng)沙運動進(jìn)行模擬的精度將更加可靠。