施文奎， 沈雁鳴， 陳堅(jiān)強(qiáng)

(中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng) 621000)

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

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

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

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

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

1 空間變光滑長(zhǎng)度SPH方法

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

Wij=W(xi-xj,h)

(1)

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

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

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

圖1 質(zhì)量初始化

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

首先在粒子搜索時(shí)采用光滑長(zhǎng)度的算術(shù)平均值來(lái)確定相互作用粒子對(duì)，即：

(2)

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

Wij=W(rij,hij),

(3)

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

(4)

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

(5)

1.2 鏈表搜索法

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

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

1.3 粒子擴(kuò)散分布模型

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

圖2 粒子擴(kuò)散分布模型

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

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

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

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

2.1 二維氣泡上浮

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

圖3 二維氣泡上浮模型

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

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

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

2.2 楔形體入水沖擊

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

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

圖6 實(shí)驗(yàn)參數(shù)描述

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

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

圖8 時(shí)間無(wú)關(guān)性驗(yàn)證Fig.8 Time independent verification

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

圖9 垂向加速度、角加速度、速度驗(yàn)證

(a) 垂向加速度

(b) 角加速度

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

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

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

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

3 結(jié) 論

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

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

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