梁嵐博，金阿芳，聞騰騰

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

風(fēng)沙運動為主要標(biāo)志的荒漠化及其帶來的風(fēng)沙環(huán)境是全球最突出的環(huán)境問題之一。風(fēng)沙流活動對農(nóng)業(yè)、牧業(yè)、沙漠公路交通、沙漠石油開采等造成了極大的危害，沙害造成的經(jīng)濟(jì)損失高達(dá)數(shù)百億元人民幣。而風(fēng)沙兩相流是一門既古老又年輕的科學(xué)，是認(rèn)識風(fēng)沙運動規(guī)律的重要基礎(chǔ)，同時也是開展沙漠化治理和促進(jìn)沙漠化區(qū)域可持續(xù)發(fā)展的科學(xué)依據(jù)，因此風(fēng)沙兩相流的研究具有重要的理論和實踐意義[1]。

風(fēng)沙兩相流的研究方法主要有風(fēng)洞模擬、野外觀測和數(shù)值模擬。風(fēng)洞試驗觀測是研究風(fēng)沙運動最常用的手段，通常采用非接觸測量儀器如粒子圖像測速儀（particle image velocimeter，PIV）、相位多普勒粒子分析儀（phase Doppler particle analyzer，PDPA）、高速攝影儀等，可以對沙粒的微細(xì)運動和風(fēng)洞輸沙的時空變化進(jìn)行分析。野外觀測更多采用風(fēng)速儀和集沙儀以及各種類型的輸沙傳感器等對風(fēng)沙運動及輸沙量進(jìn)行測量[2]。隨著計算機(jī)硬件和軟件以及計算方法的發(fā)展，數(shù)值模擬近年來在風(fēng)沙兩相流研究中的作用越來越重要。數(shù)值模擬方法便于對風(fēng)沙運動進(jìn)行定量模擬和預(yù)測，易于控制參數(shù)，相對于風(fēng)洞實驗和野外觀測有成本較低和便于實現(xiàn)的優(yōu)點。

目前，流體數(shù)值模擬研究包括基于網(wǎng)格的歐拉法和基于粒子的拉格朗日法。在基于歐拉網(wǎng)格的數(shù)值方法中如FDM，首先需要在問題域上生成網(wǎng)格，要在不規(guī)則或者很復(fù)雜的幾何形狀上構(gòu)造規(guī)律的網(wǎng)格是一件非常困難的事情，在流體模擬的整個過程中還會出現(xiàn)很多如大變形、變形邊界、自由表面和物質(zhì)交界面的問題，給基于網(wǎng)格的數(shù)值模擬方法帶來了很大的挑戰(zhàn)。拉格朗日無網(wǎng)格法主要有SPH、MPS和DEM等方法，以上粒子之間不需要進(jìn)行網(wǎng)格劃分（不是預(yù)先定義好的網(wǎng)格系統(tǒng)），是一系列任意分布的離散點，因此可以克服流體模擬中網(wǎng)格畸變和網(wǎng)格移動等引起的各種問題。目前，SPH方法被廣泛應(yīng)用于各種領(lǐng)域如潰壩、輪胎劃水、飛機(jī)開溝、血流血管和船的撞擊等[3-4]，且該方法更加穩(wěn)定和真實，但SPH方法運用在風(fēng)沙兩相流中還較為缺乏。文獻(xiàn)[5]首次將SPH 方法運用在風(fēng)沙流中，建立了基于SPH 方法的風(fēng)沙流控制方程，提出了模擬過程中一些關(guān)鍵參數(shù)的選擇，并使用虛粒子法阻止粒子穿透邊界。文獻(xiàn)[6]利用文獻(xiàn)[5]的算法，對沙粒的躍移運動進(jìn)行了模擬，并在模擬過程中引入了起風(fēng)和停風(fēng)模塊，分析了在風(fēng)的加載和卸載時，沙粒的運動軌跡。文獻(xiàn)[7]利用文獻(xiàn)[5]的算法，建立了具有坡面形態(tài)的計算模型，得到了坡面不同位置沙粒的躍移高度和水平距離，驗證了沙粒躍移的平均距離是單個沙波紋長度的1.16倍。文獻(xiàn)[8]針對歐拉-歐拉流體計算模型在求解氣沙兩項耦合流動問題中存在的缺陷，提出了基于SPH-FVM 耦合的方法模擬風(fēng)沙運動，用SPH 法離散沙粒相進(jìn)行求解，用有限體積法（finite volume method，F(xiàn)VM）求解氣相粒子，氣沙兩相通過拖曳力以及體積分?jǐn)?shù)等參數(shù)進(jìn)行耦合，模擬了自由來流風(fēng)作用下沙粒的運動過程以及沙丘緩慢向前運動的過程。由于以上模擬過程中需要將計算區(qū)域離散成數(shù)量龐大的單個粒子，其計算時間長和模擬效率低一直是該方法的瓶頸。

近年來，隨著科學(xué)技術(shù)的迅速發(fā)展，科學(xué)研究領(lǐng)域?qū)τ嬎銠C(jī)性能的要求也在不斷增長，更是對計算能力提出了極高的要求[9]。傳統(tǒng)計算機(jī)的單核處理器CPU，已經(jīng)遠(yuǎn)遠(yuǎn)不能滿足計算要求，即使是目前廣泛使用的多核處理器CPU，在龐大的計算量壓力下，也捉襟見肘[10]。為了滿足巨大計算需求，眾核架構(gòu)圖形處理器（GPU）的出現(xiàn)成為了必然。目前，由于GPU 超強(qiáng)的計算能力且適合計算密集型和易于并行的程序，以廣泛應(yīng)用到SPH 水體仿真算法中[11-12]，但應(yīng)用到風(fēng)沙運動中還比較匱乏。

因此本文引入了GPU 眾核架構(gòu)，并且通過對SPH串行算法的改進(jìn)，提出了一種SPH-GPU并行計算方法，用以在計算機(jī)上快速的模擬風(fēng)沙運動。

1 風(fēng)沙流的SPH描述和離散

1.1 SPH基本思想

SPH 方法是在1977 年由Gingold 和Monaghan 等人[13-14]提出，最初用于解決天體物理學(xué)問題，后來經(jīng)過其他學(xué)者的不斷改進(jìn)和完善，推廣到流體力學(xué)各個方面。該無網(wǎng)格方法是將連續(xù)的流體離散成一系列任意分布的SPH 粒子，每個粒子攜帶質(zhì)量、速度、密度等物理信息，通過求解SPH控制方程來描述整個系統(tǒng)的演變過程。

1.2 SPH基本方程

SPH方法的構(gòu)造按兩個關(guān)鍵步驟進(jìn)行：第一步為積分表示法，又稱場函數(shù)核近似法；第二步為粒子近似法[15]。場函數(shù)核近似法如下式：

式中，Ω是支持域，f(x)為三維坐標(biāo)向量x的函數(shù)，δ(x-x′)是狄拉克函數(shù)，性質(zhì)如下：

用光滑函數(shù)W(x-x′)取代f(x)積分表示式中的δ(x-x′)則：

式中，角括號＜＞在SPH的習(xí)慣用法中是核近似算子的標(biāo)記；W(x-x′)即光滑函數(shù)；h是光滑長度，用來定義光滑函數(shù)影響區(qū)域的范圍。

式（3）可轉(zhuǎn)化為支持域內(nèi)所有粒子疊加求和的離散化形式，稱為粒子近似法。在粒子i處的粒子近似形式可寫為：

其中，N是粒子i的支持域內(nèi)的粒子總數(shù)；xi、xj分別為i和j粒子的坐標(biāo)向量；ρj、mj為j粒子的密度和質(zhì)量。

1.3 N-S方程的SPH離散

1.3.1 連續(xù)密度法

對連續(xù)性方程進(jìn)行速度散度SPH近似：

1.3.2 動量方程的粒子近似法

對動量方程的梯度項直接應(yīng)用SPH 粒子近似法進(jìn)行變換可得：

1.4 鄰近粒子的搜索

在模擬風(fēng)沙運動的過程中，氣沙兩相粒子之間的位置是不斷變化的，所以需要在每個時間步內(nèi)都要重新計算每個粒子作用域中所包含的粒子，可以說鄰近粒子搜索法的選擇是非常重要的一步。

圖1（a）全配對搜索法是最簡單和直接的搜索方法。對于目標(biāo)粒子i（實心紅色圓點），要計算它到問題域中每一個粒子j(j=1,2,…,N,N是問題域中的粒子總數(shù)，圖中用藍(lán)色實心圓點表示）之間的距離，若距離小于粒子的支持域的半徑kh，則粒子j為粒子i的鄰近粒子（所有在紅色圓圈之內(nèi)的點）。在本文中因考慮對稱光滑長度，則粒子i和粒子j為一組相互作用的粒子對，即粒子i也在粒子j的支持域內(nèi)，那么，在搜索粒子j的鄰近粒子時就不需要再考慮粒子i，因此可以減少一半的搜索計算量。該搜索方法對計算域內(nèi)的所有粒子進(jìn)行，雖然這個方法的應(yīng)用較簡單，但顯然，這種方法的時間復(fù)雜度是O(N2)。

圖1 鄰近粒子的搜索Fig.1 Search for nearby particles

在運用鏈表搜索法時，在問題域上要鋪設(shè)一臨時網(wǎng)格，如圖1（b）所示。將每個粒子都分布在網(wǎng)格單元內(nèi)，并通過簡單的存儲規(guī)則將每個網(wǎng)格內(nèi)的所有粒子連接起來。網(wǎng)格單元的大小與粒子支持域的大小相同，此網(wǎng)格只是用來劃分區(qū)域，不參與計算。即，若粒子i的支持域大小為kh，則網(wǎng)格單元的大小也為kh。那么，對于給定的粒子i，其鄰近粒子只能在同一網(wǎng)格或相鄰網(wǎng)格里，由于本文的問題域是二維的，所以相鄰網(wǎng)格的數(shù)量為8。若每個單元內(nèi)的平均粒子數(shù)量很小，則鏈表搜索法的復(fù)雜度為O(N)。

上述算法在構(gòu)建相鄰粒子對時所使用到的輸入數(shù)據(jù)僅與當(dāng)前粒子（沙粒、氣體以及虛粒子）的空間位置有關(guān)，與其他變量無關(guān)。每次計算輸出的結(jié)果是粒子對的信息，并不改變當(dāng)前時刻的粒子位置數(shù)據(jù)。在構(gòu)建粒子對循環(huán)函數(shù)體中計算出的結(jié)果不影響輸入數(shù)據(jù)，那么函數(shù)體中每次迭代的過程都不相互依賴。因此任一時刻某個粒子在尋找相鄰粒子對時都是相互獨立的，也就是說讓GPU中的每一個線程攜帶粒子空間位置信息進(jìn)行計算，每個線程執(zhí)行時的區(qū)別僅是輸入位置數(shù)據(jù)的不同，而線程的計算指令保持一致，符合GPU大規(guī)模并行計算的要求。

2 計算條件

文中主要是基于氣沙兩相耦合之間的受力機(jī)制來模擬沙粒的運動過程。模擬過程中涉及的各類參數(shù)及運行環(huán)境分別為表1和表2所示。由于在邊界上的粒子積分時會被邊界截斷，會導(dǎo)致求解結(jié)果出錯，所以在計算區(qū)域上部和下部采用固體邊界，進(jìn)出口采用周期性邊界。在計算區(qū)域頂部和沙粒的底部分別設(shè)置一層虛粒子，從而阻止粒子非物理穿透邊界如圖2所示。初始床面上的風(fēng)速廓線公式：

圖2 計算區(qū)域初始狀態(tài)示意圖Fig.2 Schematic diagram of initial state of calculation area

表1 沙相、氣相物性和計算參數(shù)Tabel 1 Sand，gas phase physical properties and calculation parameters

表2 計算環(huán)境Tabel 2 Computing environment

式中，μz為床面高程z處的風(fēng)速；z0為沙床面粗糙度，z0=Ds/30。

3 SPH方法在CUDA計算架構(gòu)下的實現(xiàn)

為了實現(xiàn)并行計算，需要對SPH 串行程序進(jìn)行分析和修改，首先尋找程序中的熱點函數(shù)（程序中耗時最多的函數(shù)），此外還要考慮數(shù)據(jù)的依賴性（前面的計算結(jié)果緊接著被后面的函數(shù)用到）分析哪部分程序適合串行和并行計算等等。最后用C語言編寫主流程，采用CUDA C 語言調(diào)用GPU 設(shè)備實現(xiàn)并行計算，使改進(jìn)后的程序，能使CPU與GPU協(xié)同合作完成任務(wù)，最大限度的發(fā)揮計算機(jī)的并行處理能力，從而提高SPH數(shù)值計算效率。

如圖3 所示模擬過程分為4 個部分：讀取需要模擬的數(shù)據(jù)、搜尋鄰居粒子、計算粒子間作用力以及粒子位置的更新。CUDA 程序開始計算時，CPU 程序占主導(dǎo)地位，主機(jī)端創(chuàng)建輸入和輸出數(shù)組，為輸入數(shù)據(jù)和結(jié)果提供存儲空間，將已有的粒子信息文件（初始的位置、速度、密度等屬性值）讀取到CPU 內(nèi)存里；隨后CPU加載cudaMalloc函數(shù)分配GPU顯存和cudaMemcpy函數(shù)把主機(jī)內(nèi)存里的粒子信息拷貝至GPU 顯存里，接下來CPU 加載核函數(shù)給GPU 做計算，GPU 計算時CPU阻塞不執(zhí)行其他任務(wù)，核函數(shù)在設(shè)備端進(jìn)行相鄰粒子搜索，搜索之后計算粒子屬性值，在計算粒子屬性值時要考慮數(shù)據(jù)之間的依賴性，具有依賴性的數(shù)據(jù)不適合進(jìn)行并行計算，否則會導(dǎo)致計算結(jié)果不準(zhǔn)確甚至程序崩潰；最后進(jìn)行粒子位置更新，更新完成后，跳至鄰居粒子搜索進(jìn)行循環(huán)執(zhí)行，在完成所有計算任務(wù)之后，才將所得結(jié)果拷貝到內(nèi)存中，這樣就減少了內(nèi)存和顯存數(shù)據(jù)傳輸所花費的時間，提高了并行計算的效率。這個模擬過程中還要考慮合理利用設(shè)備資源，比如使SM 始終保持計算狀態(tài)，合理的分配線程，選擇合適的存儲器等。

圖3 GPU加速的SPH計算流程Fig.3 Accelerated SPH calculation process through GPU

4 CPU-GPU并行計算結(jié)果與分析

4.1 SPH串行算法熱點分析

并行算法實現(xiàn)之前，要對整個串行程序進(jìn)行熱點分析，在進(jìn)行并行加速時，要考慮串行程序中計算時間很長的程序段，也就是熱點程序。如果首先對熱點程序進(jìn)行并行實現(xiàn)，那么對整個程序計算效率的提升是非常巨大的。上結(jié)所說模擬過程中耗時的部分主要有3個：搜尋鄰居粒子、計算粒子間作用力以及粒子更新。

如圖4是模擬25 000步得到的結(jié)果，從圖中可以很明顯地看出，整個程序耗時最大的部分在于搜尋鄰居粒子上，這部分占據(jù)了整個程序一半以上的時間，是串行程序的主要熱點所在。計算粒子間作用力，花費了接近40%的時間，也是熱點程序的部分。相比較而言，粒子信息更新在總時間中所占據(jù)的比例僅僅2%。以上可得，在進(jìn)行并行算法實現(xiàn)時，應(yīng)該首先分析搜尋鄰居粒子這段程序上，尋找更好的粒子搜尋方法，其次在優(yōu)化計算粒子間相互作用力的程序，最后對于粒子更新可以選擇保留串行計算，也可以進(jìn)行并行計算。

圖4 SPH串行熱點程序Fig.4 SPH serial hotspot program

4.2 不同搜索方法的串行SPH計算時間

由上一節(jié)得出在搜尋鄰近粒子上的耗時占用了整個程序的大部分時間，因此表3比較了全配對法和鏈表法所占的耗時，由于時間復(fù)雜度由O(N2)降低到O(N)，可見鏈表法是全配對搜索法計算速度的2.3倍。

表3 不同搜索法的SPH計算時間Tabel 3 SPH calculation time for different search methods

4.3 GPU并行計算模型驗證

4.3.1 沙粒群運動的時空變化

本算例重點討論流體起動，不考慮碰撞和沖擊產(chǎn)生的瞬時沖擊力的作用。

如圖5 給出了在SPH-GPU 并行計算下時間步為1 500步圖（a），5 100步圖（b），8 400步圖（c），9 900步圖（d），19 200 步圖（e），21 600 步圖（f）的計算結(jié)果。從圖中可知，沙粒的起動是在重力、氣流拖曳力等各種力共同作用的結(jié)果，計算結(jié)果較好地描述了沙粒群的運動規(guī)律。如圖（a）、（b）所示，在氣流的作用下，平坦沙床上層沙粒產(chǎn)生加速度，當(dāng)風(fēng)速足夠大時，沙粒從地面跳起開始運動。如圖（c）、（d）中橢圓區(qū)域所示，隨著時間的推移，當(dāng)時間步進(jìn)行到8 400步時，沙床此時開始發(fā)生松動或堆積，沙床中有些沙粒會發(fā)生離開沙床起躍的趨勢，也有些沙粒其速度朝下，這些粒子會造成沙床更加緊實，沙床的松動和緊實是造成沙床產(chǎn)生沙波紋的直接原因之一。如圖（e）、（f）所示，起跳沙粒在氣流中不斷的加速，又在重力作用下有了明顯的下落。當(dāng)沙床面上起跳的沙粒數(shù)與落回床面的沙粒數(shù)相差不大時，風(fēng)沙運動進(jìn)入了比較穩(wěn)定的狀態(tài)。

圖5 沙粒群運動的時空變化Fig.5 Temporal and spatial changes of movement of sand particles

4.3.2 沙粒的運動軌跡

如圖6為床面上典型沙粒的躍移運動軌跡，從圖中可以看到，沙粒的運動軌跡具有上升段比較急劇，降落段比較平緩的特點，都呈現(xiàn)出一種非對稱的拋物線形狀，觀察到靠近右周期邊界的躍移軌跡的下落段并沒有顯示完整，這是由于進(jìn)出口采用周期性邊界，右邊界出去的粒子從左邊界進(jìn)來，以上沙粒運動軌跡與文獻(xiàn)[16]風(fēng)洞實驗中高速攝影記錄的沙粒躍移運動相符。分析這是由于沙粒在躍移過程中始終受到氣流拖曳力的作用，導(dǎo)致了在水平方向的分速度越來越大，當(dāng)沙粒飛到一定高度，又在氣流阻力和重力的作用下，迅速向下運動，因此導(dǎo)致降落段的水平位移較大。

圖6 典型沙粒躍移軌跡Fig.6 Typical sand grain trajectories

如圖7 表示沙粒運動中具有拋物線形狀但同時又有些變異的軌跡，其躍移軌跡并不是非常光滑的拋物線形狀，而是在轉(zhuǎn)折處出現(xiàn)尖角，即沙粒會急劇降落，造成此現(xiàn)象的原因可能是局部氣流縱向脈動較大，也可能是沙粒在空中與其他沙粒發(fā)生碰撞后獲得了一個向下的沖量。

圖7 變異的尖角軌跡Fig.7 Mutated sharp trajectories

4.4 GPU-CPU計算效率的對比

如圖8 顯示了粒子數(shù)分別為2 500、10 000、40 000來計算不同粒子數(shù)下GPU與CPU的效率，可以看出，使用不同的鄰近粒子搜索法在GPU中的計算效率遠(yuǎn)高于在CPU 中的計算效率。GPU 并行鏈表搜索法相較于CPU 鏈表搜索法最高可以獲得3.1 倍的加速比，而相對于CPU全配對搜索法最高可以獲得20倍的加速比。當(dāng)粒子數(shù)目為2 500 時，GPU 并行計算的加速效果并不明顯，其原因是GPU的單線程計算能力遠(yuǎn)弱于CPU，隨著粒子數(shù)目的增加GPU多線程并行計算能力的優(yōu)勢就發(fā)揮了出來，能夠帶來更大的加速比。

圖8 CPU-GPU計算時間Fig.8 CPU-GPU calculation time

為了更好地理解CPU 與GPU 的加速機(jī)制，以10 000 粒子為例，統(tǒng)計了搜尋鄰居粒子、計算粒子間作用力以及粒子更新在計算時占總時間的比例。如圖9所示，GPU并行計算相較于CPU計算時，搜尋鄰居粒子和粒子更新所占時間比例都減少了，然而計算粒子間作用的時間比例增加了。綜合以上，可以推測出搜尋鄰居粒子和粒子更新相對于計算粒子間作用力有更好的并行計算效率。

圖9 SPH各部分占總時間的比例Fig.9 Proportion of each parts of SPH in total time

5 結(jié)論

本文首先對串行SPH程序進(jìn)行分析和修改，最后采用SPH-GPU大規(guī)模并行加速技術(shù)對風(fēng)沙流動過程進(jìn)行數(shù)值模擬，結(jié)果表明：

（1）對SPH 串行算法代碼進(jìn)行了熱點程序分析，得到在搜尋鄰居粒子和計算粒子間作用力時耗時最長，因此優(yōu)先優(yōu)化這兩部分程序段。由于搜尋鄰居粒子是最耗時的部分，比較了全配對搜索法和鏈表搜索法的計算效率，得到鏈表法是全配對搜索法計算速度的2.3 倍。

（2）對SPH-GPU 并行計算模型進(jìn)行驗證。宏觀上得到了沙粒群在各種力的作用下有起跳、上升和回落三個階段的時空變化規(guī)律以及沙波紋形成過程，微觀上得到了典型的沙粒躍移軌跡和變異的尖角軌跡。

（3）比較了CPU 與GPU 的計算效率，得到GPU 大規(guī)模并行計算的執(zhí)行時間遠(yuǎn)遠(yuǎn)低于CPU 的執(zhí)行時間。在粒子數(shù)較少時GPU 加速效果并不大，隨著粒子數(shù)的增加，加速效果越加明顯，最高可以獲得20 倍的加速比。為了更好地理解CPU 與GPU 的加速機(jī)制，比較了搜尋鄰居粒子、計算粒子間作用力以及粒子更新在計算時占總時間的比例，推測出搜尋鄰居粒子和粒子更新相對于計算粒子間作用力有更好的并行計算效率。