施文奎， 沈雁鳴， 陳堅強

(中國空氣動力研究與發(fā)展中心，四川 綿陽 621000)

SPH方法[1-2]是一種無網(wǎng)格拉格朗日數(shù)值方法，在求解模擬多介質(zhì)、多相、非定常、流固耦合、界面變形和強非線性等問題上具有獨特優(yōu)勢。此方法最早用于天體物理領(lǐng)域中，隨著方法不斷發(fā)展和改進，其在流體和固體力學(xué)領(lǐng)域應(yīng)用也越來越廣。

在SPH方法中，光滑長度h是非常重要的一個參數(shù)，是核函數(shù)的重要變量之一，直接影響到求解工程問題時的計算效率和精度。傳統(tǒng)SPH方法大多采用恒定一致光滑長度算法[3]，這就限制了整個流場區(qū)域的離散粒子為均勻分布模式。伴隨著求解問題越來越復(fù)雜，其計算域的粒子數(shù)需求也越來越大，這導(dǎo)致了計算效率嚴(yán)重下降，計算成本過于巨大。

為提高SPH方法的計算效率及空間分辨能力，需要根據(jù)流場區(qū)域重要性不同而有針對性地布置粒子分布密度。類似于有限差分法中的網(wǎng)格分布方式，在重點關(guān)心區(qū)域如沖擊位置分布較密粒子，而遠(yuǎn)離沖擊區(qū)域分布相對稀疏粒子，并給每個粒子配置獨立的光滑長度[4]和質(zhì)量[5]。

國內(nèi)這方面的相關(guān)研究相對較少，大多集中在應(yīng)用層面。強洪夫等[12]提出了完全變光滑長度SPH法，并成功應(yīng)用于高能炸藥爆轟過程的模擬。這里提出的方法和文獻[8-9]是相似的，區(qū)別在于更直接地考慮了變光滑長度的影響以及采用了對稱形式的近似核函數(shù)。此外蘇文周[13]、強洪夫等[14]也基于上述完全變光滑長度方法求解了不同類型的侵徹問題。

以往關(guān)于變光滑長度SPH方法的研究大多是為了提高計算精度，因而將粒子光滑長度與密度進行了關(guān)聯(lián)，并隨著時間推進而不斷改變，且其應(yīng)用領(lǐng)域也大多集中在天體物理、激波管和爆轟等領(lǐng)域。本文的研究目的是通過設(shè)定合適的粒子初始分布，相應(yīng)地賦予每個粒子獨立的光滑長度和質(zhì)量，并提出適合變光滑長度的鏈表搜索法，以此提高SPH方法的計算效率和空間分辨力，并將該方法成功應(yīng)用于模擬氣液、液固耦合等典型問題的求解。此方法也可為求解三維入水沖擊等復(fù)雜工程問題提供技術(shù)支撐。

1 空間變光滑長度SPH方法

SPH方法中，對于給定的粒子i，通過應(yīng)用粒子近似法，粒子i處的函數(shù)及其導(dǎo)數(shù)的表達(dá)式為

Wij=W(xi-xj,h)

(1)

式中，Wij表示相互作用粒子對的光滑核函數(shù)，h是定義光滑函數(shù)W影響區(qū)域的光滑長度，直接影響到計算精度和效率。光滑函數(shù)的支持域為x-x′ ≤κh，κ是與點x處光滑函數(shù)相關(guān)的常數(shù)。

1.1 粒子初始參數(shù)賦值及作用對稱化

當(dāng)根據(jù)計算域的重要程度確定粒子初始分布后，同種類型粒子的密度保持一致，水粒子壓力考慮靜水平衡后確定，而質(zhì)量則和粒子所占面積或體積成正比[15]，如圖1所示。此外粒子間距Δx是非均勻的，故將每個粒子獨立的光滑長度和其自身的Δx比值取為定值，從而保證每個粒子支持域內(nèi)粒子數(shù)量基本不變。此外,為保證支持域內(nèi)的粒子數(shù)量適當(dāng)，選取光滑長度h=1.23Δx[16]。

圖1 質(zhì)量初始化

當(dāng)粒子非均勻分布時，粒子i和粒子j光滑長度可能不相等。此時粒子i的影響域可能包含粒子j，但粒子j的影響域則不一定包含粒子i，這不符合牛頓第三定律。因此需要采取特定的處理以保證粒子作用的對稱性。

首先在粒子搜索時采用光滑長度的算術(shù)平均值來確定相互作用粒子對，即：

(2)

在確定了相互作用粒子對后，將求得的對稱光滑長度代入核函數(shù)中，即可求得對應(yīng)的核函數(shù)及其導(dǎo)數(shù)

Wij=W(rij,hij),

(3)

將求得的核函數(shù)及其導(dǎo)數(shù)代入SPH方程組中即可得到相應(yīng)的空間變光滑長度SPH方程組

(4)

SPH方法中很重要的一個性質(zhì)是將所有流體都考慮為微可壓的，即采用如下形式的狀態(tài)方程[17]

(5)

1.2 鏈表搜索法

傳統(tǒng)鏈表搜索法是在問題域上鋪設(shè)一臨時均勻網(wǎng)格，網(wǎng)格單元空間大小與支持域大小一致。此時對于給定粒子i而言，其相鄰粒子只能在同一網(wǎng)格單元內(nèi)或者在緊密相鄰的單元內(nèi)，算法的復(fù)雜度階數(shù)為O(N)。但當(dāng)粒子的光滑長度不一致時，如網(wǎng)格單元空間大小與特定粒子i關(guān)聯(lián)的話，則會發(fā)生漏掉粒子對的情況。此時鏈表搜索法必須加入變光滑長度處理。對于密度變化不劇烈如兩相流、入水沖擊等問題，光滑長度不隨時間改變，此時只需將網(wǎng)格單元空間大小與初始最大光滑長度粒子的支持域大小保持一致即可。對于密度變化劇烈如爆炸等問題，光滑長度隨時間改變，此時則需要在劃分網(wǎng)格前將粒子遍歷一遍找到最大光滑長度粒子。

在相同粒子數(shù)的情況下，加入變光滑長度處理后的鏈表搜索法不可避免地會增加搜索時間，但其效率仍比直接搜索法要高很多，且此時其算法的復(fù)雜度階數(shù)仍為O(N)。

1.3 粒子擴散分布模型

傳統(tǒng)SPH方法大都采用區(qū)域一致恒定光滑長度算法，粒子為均勻分布方式。本文采用的是粒子擴散分布模型，在重點關(guān)心區(qū)域布置較密粒子，而在遠(yuǎn)離重點區(qū)域所布置的粒子逐漸稀疏。氣泡上浮算例和楔形體入水沖擊算例粒子分布示意圖,見圖2。

圖2 粒子擴散分布模型

氣泡上浮算例中，氣泡直徑為2R，水域大小為10R×6R。此算例中氣泡垂直上升，為左右對稱算例。氣泡上升階段的形狀、位置是我們需要重點關(guān)注的信息。故設(shè)計粒子初始分布時將計算域分為三部分，中間段(寬為2R)粒子采用傳統(tǒng)均勻分布方式，粒子間距Δx/R=0.05。兩側(cè)粒子間距從中間開始逐步增大，最大處橫向粒子間距為Δx/R=0.144。此時相對全均勻分布模型粒子數(shù)減少了約1/4。

非對稱楔形體自由入水沖擊算例中，沖擊載荷峰值、角加速度峰值等是我們重點關(guān)心的信息。因此設(shè)計粒子初始分布時將水域分為上下兩部分(上∶下=1∶2)，上層粒子采用均勻分布，粒子間距為0.005。下層粒子從分界處開始逐步增大，最大處縱向粒子間距為0.015。此時相對全均勻分布模型粒子數(shù)減少了約1/4。

2 計算結(jié)果及分析

計算中壓力P通過采用了人工壓縮率的狀態(tài)方程(式(5)求解，水的聲速值取14(gR)0.5。核函數(shù)采用Gussian核函數(shù)。時間推進采用Leap-Frog格式，邊界處理采用Liu等[18]提出的第一類粒子模型。下面通過氣泡上浮和楔形體入水沖擊算例，主要展示粒子擴散分布模型對計算效率的改進效果，并分析了空間變光滑長度算法對由粒子非均勻分布引起的不穩(wěn)定性的抑制作用。

2.1 二維氣泡上浮

氣泡上浮是一種典型的氣液兩相自由表面流動。圓形氣泡在浮力作用下逐漸上浮，伴隨著氣泡變形甚至破裂。過程中壓差阻力、表面張力作用較為復(fù)雜。選取的計算模型見圖3，初始時刻，氣泡半徑R為0.3 m的氣泡，位于寬為6R，高為10R的水域中。水與空氣密度比為：ρX∶ρY=1 000∶1。

圖3 二維氣泡上浮模型

首先采用傳統(tǒng)SPH方法對均勻粒子分布模型進行模擬，此時粒子間距均為Δx=0.05R，光滑長度為h=1.23Δx，粒子總數(shù)為24 940個。時間步長取5×10-6s，采用4核并行計算。同時為保證交界面處的壓力匹配，狀態(tài)方程中γX=7，γY=1.4，且取水的聲速值為cX/(gR)1/2=14，空氣的聲速值為cY/(gR)1/2=198。XSPH方法中的參數(shù)ε取為0.1，表面張力項采用文獻[19]中給出的形式。圖 4給出了計算結(jié)果與Level-Set計算結(jié)果[20]的比較?？梢钥闯?，氣泡形狀和位置在不同時刻都吻合得較好，驗證了程序模擬氣液兩相流動的有效性。

在此基礎(chǔ)上，為提高計算效率，采用粒子擴散分布模型(圖 2(a))。此時粒子數(shù)為19 120個，粒子數(shù)減少了約1/4。其它計算參數(shù)選擇和上面一致。圖5展示了采用粒子擴散分布模型時空間變光滑長度算法對結(jié)果的改進效果。三角形的點是圖4中傳統(tǒng)SPH方法計算得到的結(jié)果。圓圈的點是利用粒子擴散分布模型恒定光滑長度算法計算得到的結(jié)果，效果很差。正方形的點是擴散分布模型空間變光滑長度算法計算得到的結(jié)果，可以看出氣泡位置和形態(tài)都與粒子均勻分布模型計算得到的結(jié)果吻合得很好。

表1給出了粒子均勻分布模型和粒子擴散分布模型計算效率的比較。為避免多核并行帶來的干擾，這里采用單核進行測試。此外由于不同粒子搜索方法對于效率的對比會產(chǎn)生很大的偏差，因此這里的計算時間不包括粒子搜索時間。結(jié)果表明，計算時間節(jié)約1/4左右，大幅提高了計算效率。

2.2 楔形體入水沖擊

楔形體入水沖擊是一種典型的流固耦合問題。帶初始傾斜角的入水沖擊問題耦合了大變形、多自由度等現(xiàn)象，其模擬難度更大。選取的計算模型及相關(guān)參數(shù)見圖6，更詳細(xì)的楔形體實驗參數(shù)可參考文獻[21]。本算例中，水域大小設(shè)為長×高=3.6 m×1.8 m。

表1 不同粒子分布模型計算時間比較(不包含搜索時間)

圖6 實驗參數(shù)描述

首先采用傳統(tǒng)SPH方法對粒子均勻分布模型進行模擬，此時粒子間距均為0.005，光滑長度均為h=1.23Δx，粒子總數(shù)為265 277個。時間步長取10-5s，采用4核并行計算。狀態(tài)方程中水的聲速值取14(gH)0.5,H為水深1.8 m。固壁邊界壓力采用SPH粒子近似法得到，流固耦合模塊計算采用文獻[22]中的方法。在上述參數(shù)基礎(chǔ)上，首先進行了粒子無關(guān)性(圖7)和時間無關(guān)性驗證(圖8)，結(jié)果表明粒子間距和時間步長的選取是有效的。圖9給出了垂向加速度、角加速度和楔形體速度計算結(jié)果與實驗數(shù)據(jù)對比圖。由圖可得，程序?qū)α鞴恬詈虾鸵好孀冃蔚哪M是有效的。

圖7 粒子無關(guān)性驗證Fig.7 Particle independence verification

圖8 時間無關(guān)性驗證Fig.8 Time independent verification

在此基礎(chǔ)上，為提高效率，采用粒子擴散分布模型，見圖 2(b)。此時粒子數(shù)為199 029個，粒子數(shù)減少了約1/4，時間步長仍取10-5s。圖10給出了采用粒子擴散分布模型時恒定光滑長度算法和空間變光滑長度算法的結(jié)果對比。由圖可得恒定光滑長度算法得到的計算結(jié)果較差，這是由于粒子非均勻分布導(dǎo)致強不穩(wěn)定性造成的。但當(dāng)采用了空間變光滑長度技術(shù)后，首先其加速度和角加速度與均勻分布粒子模型計算結(jié)果吻合得都很好；其次其流場紊亂程度降低降慢(圖11)。這一方面是給每個粒子單獨配置光滑長度后，其支持域內(nèi)的粒子數(shù)量能夠保證基本一致；另一方面則是實現(xiàn)了粒子相互作用對稱化，保證了牛頓第三定律。

圖9 垂向加速度、角加速度、速度驗證

(a) 垂向加速度

(b) 角加速度

(a) t=0.377 s(左、中和右分別為均勻模型，非均勻變h，非均勻恒定h)

(b) t=0.407 s(左、中和右為均勻模型，非均勻變h，非均勻恒定h)

表2則給出了粒子均勻分布模型和粒子擴散分布模型計算效率的比較。和氣泡上浮算例類似，這里采用單核測試，計算時間不包括粒子搜索時間。結(jié)果表明計算時間節(jié)約1/4左右，大幅提高了計算效率。

表2 不同粒子分布模型計算時間比較(不包含搜索時間)

3 結(jié) 論

為提高計算效率，類似于有限差分法中的網(wǎng)格分布方法，根據(jù)計算域中物理量的變化劇烈程度匹配相應(yīng)疏密的粒子分布，構(gòu)建了粒子擴散分布模型。然后給每個粒子單獨配置相對應(yīng)的光滑長度，通過對粒子對的光滑長度取算術(shù)平均值來實現(xiàn)粒子作用對稱化，并對鏈表搜索法添加了變光滑長度處理。通過氣泡上浮和楔形體入水沖擊兩個典型算例驗證了粒子擴散分布模型及空間變光滑長度算法的有效性，結(jié)論表明：

(1) 采用粒子擴散分布模型能夠有效減少粒子數(shù)量，但會引起強不穩(wěn)定性，導(dǎo)致結(jié)果出錯。而空間變光滑長度算法可以很好地抑制上述由于粒子分布不均勻?qū)е碌膹姴环€(wěn)定性，在重點關(guān)心的氣泡信息和楔形體沖擊要素等方面都取得很好的模擬效果。

(2) 合理粒子分布方式下，利用空間變光滑長度算法可以在保證結(jié)果準(zhǔn)確性的同時大大降低計算成本，為三維復(fù)雜工程問題的求解打下基礎(chǔ)。但在粒子分布方式的優(yōu)化、標(biāo)準(zhǔn)化及抑制變光滑長度引起的不穩(wěn)定性方面需要進一步深入研究。