許曉陽，王斯棋

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

1 引言

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

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

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

2 數(shù)學模型與方法

2.1 控制方程

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

(1)

(2)

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

(3)

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

(4)

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

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

2.2.1 控制方程的離散

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

(5)

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

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

(6)

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

(7)

對函數(shù)空間導數(shù)進行適當?shù)臄?shù)學處理，利用式(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用于防止粒子相互靠近時產(chǎn)生的數(shù)值振蕩。

2.2.2 邊界處理

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

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

(10)

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

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

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

2.2.3 時間積分

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

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

3.1 算法流程

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

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

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

3.2 評價并行算法的參數(shù)

加速比和并行效率是衡量并行算法性能的2個關(guān)鍵參數(shù)。

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

4 并行計算模型驗證

4.1 有效性驗證

為了驗證基于粒子分解的SPH并行算法模擬自由表面流問題的有效性，對二維潰壩問題進行數(shù)值模擬。圖3給出了二維潰壩初始狀態(tài)模型示意圖，其幾何尺寸與文獻[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，所用流體粒子數(shù)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給出了該算例使用不同進程數(shù)運行4萬個時間步的并行結(jié)果。從表1可以看出，當使用2個或4個進程時，可以獲得較好的加速比和并行效率，但當所使用的進程數(shù)目增加到8個時，加速比和并行效率出現(xiàn)了一定程度的下降。這是由于模擬所用粒子數(shù)較少，當所用進程數(shù)增加時，分配給每個進程的粒子數(shù)相應(yīng)減少，這造成需要通信的數(shù)據(jù)占總數(shù)據(jù)的比例增大，通信量增加，因此加速比和效率出現(xiàn)了下降。

Table 1 Parallel performance analysis of dam-break flow

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

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

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。粒子之間初始間距設(shè)置為δ0=0.000105 m，所用粒子總數(shù)設(shè)置為N=1071521，其中流體粒子總數(shù)設(shè)置為825 421，邊界粒子總數(shù)設(shè)置為58 801，固壁外虛粒子總數(shù)設(shè)置為187 299。時間步長Δt=5.0×10-7s，以保證數(shù)值穩(wěn)定性。

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

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

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

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

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

5 結(jié)束語

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

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

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

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

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