基于粒子分解的SPH并行算法研究與應用*

許曉陽，王斯棋

(西安科技大學計算機科學與技術學院，陜西 西安 710054)

1 引言

自由表面流現象廣泛存在于自然界和工業生產中，如注塑成型、水利工程等。處于流動過程的自由面形狀復雜，且在演變中可能產生水花四濺、靠近邊界處水流猛烈變形、反彈以及與下方水體交融等多種復雜物理現象，是一個強非線性復雜問題。因此，高效、準確地模擬這一流動過程具有重要的理論價值和實用意義?；趥鹘y網格的數值方法如有限差分、有限元法等在解決此類強非線性自由表面流問題時，需要運用額外界面追蹤技術，實施過程復雜。

光滑粒子流體動力學SPH(Smoothed Particle Hydrodynamics)方法是一種拉格朗日型無網格方法，適用于模擬自由表面流動。此方法是由Lucy[1]和Gingold等人[2]在1977年首次提出的。與基于網格的數值方法相比，SPH方法完全獨立于網格，且具有Lagrangian特性和質點特性、自適應特性等優點，因此非常適合大變形、自由表面流動等復雜界面問題的數值模擬。1994年，Monaghan[3]首次將SPH方法應用于自由面的不可壓縮建模，基于可壓縮假設構建了一個計算簡單的不可壓縮流體模型。近年來，SPH方法又被成功用于求解不可壓縮流[4,5]、多相流[6,7]、傳熱[8,9]和粘彈性流[10,11]等問題。值得注意的是，相比于網格類方法，SPH方法在計算過程中計算量大、耗時長，因此有必要對SPH程序進行并行化處理。目前，王裴等人[12]采用固定空間區域方法實現了三維微噴射和斜侵徹的并行SPH模擬。Cherfils等人[13]設計了基于區域分解的SPH并行算法，并對二維水柱倒塌過程中的并行性能進行了研究與分析。吳建松等人[14]借助GPU并行加速技術，應用SPH方法對復雜階梯流問題進行了數值模擬。范小康等人[15]利用基于可伸縮矢量擴展SVE(Scalable Vector Extension)的單指令多數據SIMD(Single Instruction Multiple Data)結構向量優化方法對SPH數量級并行進行了探索，獲得了明顯加速效果。梁嵐博等人[16]在CUDA軟硬件平臺上，建立了SPH-GPU并行加速二維氣沙兩相耦合模型，結果表明該方法能夠進一步應用在風沙流數值模擬中。

本文基于消息傳遞接口MPI(Message Passing Interface)并行程序設計平臺，以C++語言作為算法實現的編程語言，設計了基于粒子分解的SPH并行算法。該算法將所有粒子平均分配到各個進程進行計算，每個時間步通信僅調用一次發送、接收和廣播函數，因此易于實現且可擴展性較好。應用該并行算法對二維潰壩流和三維液滴沖擊液膜問題進行了模擬，分析了進程數、粒子數與并行效率、加速比之間的關系。當粒子數大于百萬時，最大加速比可達30以上，為進行三維大規模問題的數值模擬提供了一種高效的計算工具。

2 數學模型與方法

2.1 控制方程

在Lagrangian坐標系下，三維等溫、牛頓黏性流體的控制方程組如式(1)和式(2)所示:

(1)

(2)

其中，ρ、v和p分別表示流體密度、速度和壓力，t表示時間，η表示流體粘度，g表示重力，D/Dt表示物質導數，其定義如式(3)所示:

(3)

為封閉控制方程式(1)和式(2)，通常將不可壓縮流體視為弱可壓縮流體，即用狀態方程將流體密度的變化范圍控制在1%以內，以保證流體流動行為完全接近不可壓縮流。本文使用的狀態方程如式(4)所示[3]：

(4)

其中，γ=7是一個常數；ρ0表示參照密度，其值為1 000 kg/m3；c表示聲速，通常取為流體最大速度的10倍。

2.2 光滑粒子流體動力學(SPH)方法

2.2.1 控制方程的離散

在SPH方法中，函數f(x)在計算域Ω內的積分表達式可寫為式(5)所示:

(5)

其中，x表示位置矢量，W(·)表示核函數，h表示核函數影響域的光滑長度。本文核函數采用分段三次樣條函數，此時影響域半徑為2h。

對計算域Ω內第i個粒子的位置矢量，將式(5)轉化為核函數支持域內粒子疊加求和的離散化形式,如式(6)所示:

(6)

其中，N表示粒子總數,mj表示第j個粒子的質量；Wij=W(|xij|,h);xij=xi-xj，|xij|表示第i個粒子和第j個粒子之間的距離。同理求得，函數空間導數在粒子i處的粒子近似式如式(7)所示:

(7)

對函數空間導數進行適當的數學處理，利用式(7)可推導出不同的粒子近似式，進而用于控制方程式(1)和式(2)的離散。本文選用的離散形式如式(8)和式(9)[17]所示：

(8)

(9)

其中，vij=vi-vj;ηi和ηj分別表示第i個粒子和第j個粒子的粘度;pi和pj分別表示第i個粒子和第j個粒子的壓力;ρi和ρj分別表示第i個粒子和第j個粒子的流體密度;φ=0.01h用于防止粒子相互靠近時產生的數值振蕩。

2.2.2 邊界處理

邊界處理直接影響模擬的效率和穩定性，對于SPH計算至關重要。本文在前期工作[18,19]基礎上，提出一種由壁面粒子和邊界外虛粒子組成的加強型邊界處理方法。

首先，在固壁邊界上布置一層壁面粒子，且粒子間距與流體粒子初始間距δ0相等。與Monaghan[3]的邊界方法不同，本文壁面粒子不通過施加排斥力以防御流體粒子穿透固壁。與文獻[20]和文獻[21]的方法類似，壁面粒子參與到流體控制方程的求解中。在計算過程中，壁面粒子的密度和位置不發生變化，壓力通過其支持域內流體粒子壓力的正則化插值計算得到，如式(10)所示：

(10)

其中，i表示壁面粒子，j表示與壁面粒子i相鄰的流體粒子。

其次，在固壁邊界外布置幾層虛粒子，以彌補壁面粒子的不足。虛粒子與流體粒子初始間距δ0相等，密度和位置在計算過程中保持不變。但與文獻[20]和文獻[21]的方法不同，虛粒子的速度和壓力不再通過構造流體內部偽粒子進行插值計算得到。本文中，每個虛粒子均有唯一的壁面粒子與之相連接。圖1展示了虛粒子與壁面粒子的連接關系。為符合無滑移邊界條件，壁面粒子和虛粒子的速度均設置為零。虛粒子壓力設置與相連接固壁粒子的壓力相同。

相比于文獻[20]和文獻[21]的方法，本文的邊界處理將不再需要構造流體區域靠近固壁邊界處的偽粒子，從而避免了應用式(10)對偽粒子進行正則化插值的繁瑣操作，因此可以縮短三維模擬的計算時間。

2.2.3 時間積分

由于蛙跳法具有二階精度，且對于三維問題的存儲需求量小、計算效率高，因此本文選用該方法對SPH離散方程式(8)和式(9)進行時間積分。關于蛙跳法的時間步推進公式可參閱文獻[17]。

3 基于粒子分解的SPH并行算法

3.1 算法流程

本文基于MPI并行程序設計平臺，以C++語言作為算法實現的編程語言，設計了一套基于粒子分解的SPH并行算法。

該并行算法的基本思想是把相鄰粒子間相互作用力的計算，根據各處理器計算能力各自分配一定數量的粒子，進行每一時間步的通信和并行計算。首先，輸入初始粒子信息和計算所需的其它數據，將所有粒子平均分配到各個進程：設總粒子數為N，總進程數為P，進程數標記為Z(0≤Z≤P-1)，先計算N/P和N%P，若Z>N%P，則分配給進程Z的粒子起止編號分別為Z×(N/P)+N%P和(Z+1)×(N/P)+N%P-1；否則，分配給進程Z的粒子起止編號分別為Z×(N/P+1)和(Z+1)×(N/P)+Z。其次，對流體控制方程進行并行求解。計算過程中，每一時間步通信僅調用一次發送、接收和廣播函數，因此編程易于實現，且可擴展性較好。該并行算法的另一特點是每個進程在每一時間步內均負責維護固定的某一部分粒子，并不考慮粒子實際所處物理空間，因此各進程間的負載平衡易于保證。數值算例表明，該并行算法可顯著提升SPH方法模擬三維復雜流動問題的計算能力，對其他如耗散粒子動力學、分子動力學程序并行也可提供有價值的參考。

圖2展示了該并行算法的流程。

3.2 評價并行算法的參數

加速比和并行效率是衡量并行算法性能的2個關鍵參數。

加速比Sn是指同一任務串行運行時間T1與并行運行時間Tn之比，其中n表示所用總進程數。并行效率En是指并行加速比與總進程數之比。En一般小于或等于1，越接近1說明并行加速效率越高。

4 并行計算模型驗證

4.1 有效性驗證

為了驗證基于粒子分解的SPH并行算法模擬自由表面流問題的有效性，對二維潰壩問題進行數值模擬。圖3給出了二維潰壩初始狀態模型示意圖，其幾何尺寸與文獻[13]的保持一致，即潰壩水體高度H=1 m,長度L=2 m，水槽寬度d=5.366 m,m為水槽高度。流體密度ρ=1 000 kg·m-3，粘度η=10-3Pa·s，重力加速度g=9.81 m·s-2，所用流體粒子數NF=5 000，時間步長Δt=1.0×10-4s。

Figure 3 Schematic diagram of 2D dam-break model 圖3 二維潰壩模型示意圖

Figure 4 SPH particle distribution of dam-break flow at two different times圖4 潰壩流在2個不同時刻的粒子分布圖

Figure 5 Front position of dam-break flowtime history varying varying with time圖5 潰壩流前沿位置隨時間的變化圖

表1給出了該算例使用不同進程數運行4萬個時間步的并行結果。從表1可以看出，當使用2個或4個進程時，可以獲得較好的加速比和并行效率，但當所使用的進程數目增加到8個時，加速比和并行效率出現了一定程度的下降。這是由于模擬所用粒子數較少，當所用進程數增加時，分配給每個進程的粒子數相應減少，這造成需要通信的數據占總數據的比例增大，通信量增加，因此加速比和效率出現了下降。

Table 1 Parallel performance analysis of dam-break flow

另外，為進一步分析計算規模增大時的并行性能，圖6分別展示了潰壩流在粒子數NF=5000和NF=20000 時的并行效率?？梢钥闯?，當所用流體粒子數增加到NF=20000時，使用8個進程時的并行效率較NF=5000時有明顯提高，且在85%以上。因此依據等效率可擴展性度量法[22],本文SPH并行算法具有良好的擴展性。

Figure 6 Parallel efficiency of dam-break flow at different computation scales圖6 潰壩流在不同計算規模時的并行效率

4.2 三維液滴沖擊液膜問題

接下來，本文對三維液滴沖擊液膜問題進行SPH并行模擬，其計算模型如圖7所示。其中，液滴和液膜采用相同的液體，密度ρ=1200 kg·m-3，粘度η=0.022 Pa·s。液滴直徑D=0.0042 m，液滴沖擊速度V=5.09 m·s-1。液膜長度和寬度分別為Lx=Ly=5D，液膜厚度H′=0.5D。粒子之間初始間距設置為δ0=0.000105 m，所用粒子總數設置為N=1071521，其中流體粒子總數設置為825 421，邊界粒子總數設置為58 801，固壁外虛粒子總數設置為187 299。時間步長Δt=5.0×10-7s，以保證數值穩定性。

Figure 7 Calculation model of 3D droplet impacting liquid film圖7 三維液滴沖擊液膜的計算模型

圖8給出了三維液滴沖擊液膜問題在4個不同時刻的SPH結果?？梢钥闯?，在t=0.25 ms時，液滴以一定初始速度沖擊靜止的液膜，液滴部分粒子和受沖擊液膜部分粒子相互融合，躍出了液膜表面，形成了薄片射流。在t=0.50 ms時，更多液滴粒子持續沖擊液膜，液膜沿著固體壁面向四周逐漸擴展和運動，而薄片射流也持續向上運動，形成了較明顯的“皇冠”狀水花。在t=1.00 ms時，由于慣性力的作用液滴繼續向下運動，最終完全與液膜融合在一起，而部分從“皇冠”狀水花邊緣脫落的粒子形成小液滴，最終產生水花飛濺現象。很明顯，本文SPH并行算法能夠形象逼真地捕捉液滴沖擊液膜發生的“皇冠”狀水花、水花飛濺等多種復雜物理變化。

Figure 8 SPH simulation of 3D droplet impacting liquid film圖8 三維液滴沖擊液膜的SPH模擬

表2進一步比較了在不同進程數情況下，液滴沖擊液膜問題運行一個時間步長消耗的計算時間和并行加速比?？梢钥闯觯瑢τ诹Ｗ訑荡笥诎偃f的本算例，使用串行程序(進程數為1)所消耗的運行時間為180.9 s，而使用64個進程時計算時間縮減到5.82 s，此時最大加速比可達30以上。這說明，本文基于粒子分解的SPH并行算法能顯著減少模擬所用時間，有利于進行三維大規模計算問題的數值模擬。

Table 2 Parallel performance analysis of 3D droplet impacting liquid film

5 結束語

為解決SPH方法計算量大、耗時長的問題，本文提出了基于粒子分解的SPH并行算法。通過數值模擬二維潰壩流、三維液滴沖擊液膜問題，所得結論如下所示：

(1)二維潰壩流問題的數值結果與文獻結果相吻合，驗證了本文SPH并行算法模擬自由表面流問題的有效性；

(2)當粒子數較少、進程數較多時，通信量增加，導致加速比和并行效率出現一定程度的下降。

(3)對粒子數大于百萬的三維復雜流動問題，最大加速比可達30以上。

為進一步提升SPH方法模擬復雜流動問題的計算能力，后續將開展基于GPU的SPH并行算法研究。