粒子法中自由表面粒子識別

， ， ，

(清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，教育部先進(jìn)反應(yīng)堆工程與安全重點(diǎn)實(shí)驗(yàn)室，北京 100084)

1 引 言

移動粒子半隱式法MPS(Moving Particle Semi-implicit Method)是一種基于拉格朗日Lagrange粒子的無網(wǎng)格方法，最早由Koshizuka等[1-3]提出，并隨后運(yùn)用于計(jì)算核能領(lǐng)域的一些熱工水力難題，如蒸汽爆炸、兩相流和液滴破裂等大變形問題及復(fù)雜的流固耦合問題。上述問題如果采用傳統(tǒng)的網(wǎng)格方法進(jìn)行模擬，在計(jì)算過程中容易出現(xiàn)網(wǎng)格變形擠壓，從而造成計(jì)算不收斂。而MPS方法完全克服了上述困難，采用完全的Lagrange方法進(jìn)行模擬計(jì)算，粒子之間通過核函數(shù)產(chǎn)生相互聯(lián)系，采用梯度模型和拉普拉斯模型等離散控制方程，由于粒子之間沒有拓?fù)潢P(guān)系，所以對自由表面變化非常劇烈的物理現(xiàn)象能夠進(jìn)行較好模擬，也能夠取得較好的計(jì)算結(jié)果。

由于MPS方法發(fā)展時(shí)間較短，仍存在一些缺陷。比如在傳統(tǒng)的MPS方法中，一直沒有解決如何準(zhǔn)確判定自由表面粒子的問題，而這一問題也是導(dǎo)致在MPS方法計(jì)算中壓力計(jì)算出現(xiàn)震蕩的原因。Lee等[4]發(fā)現(xiàn)，MPS方法在判定自由表面粒子方面存在缺陷，比如容易將流體內(nèi)部較稀疏區(qū)域的粒子識別成自由表面粒子；在計(jì)算過程中，流體內(nèi)部壓力變化不均勻，在局部會出現(xiàn)壓力突變，從而導(dǎo)致流體內(nèi)部壓力振蕩和由于壓力振蕩導(dǎo)致的流體表面粒子飛濺。

針對基于Lagrange方法的粒子法中自由表面粒子判定問題，其主要方法可以分為幾何法和體積法兩類。體積法方面，Randles等[5]對存在破碎的物質(zhì)運(yùn)用SPH方法進(jìn)行模擬時(shí)，通過采用重正化的核函數(shù)來改進(jìn)離散誤差，從而解決了在自由表面接觸處密度不連續(xù)的問題，改善了自由表面的判定；Lee等[6]在對比WCSPH方法和ISPH方法時(shí)，對自由表面邊界條件進(jìn)行了改進(jìn)，通過粒子的位置散度大小來確定粒子是否為自由表面粒子，該方法能判斷出大部分自由表面粒子，但還有小部分自由表面粒子無法確定；Tanaka等[7]在研究MPS方法中的壓力振蕩問題時(shí)，通過重新構(gòu)造自由表面處的粒子數(shù)密度和初始粒子數(shù)密度的方法，改進(jìn)了自由表面粒子的判別，一定程度改善了MPS計(jì)算中的壓力振蕩問題。幾何法方面，Dilts[8]在開發(fā)全新的MLSPH方法時(shí)，針對該方法在收斂性方面存在的問題，開發(fā)了全新的基于幾何方法的自由表面粒子判定方法，并通過對球形與平面碰撞的模擬驗(yàn)證了自由表面判斷的改進(jìn)；Haque等[9]借鑒Dilts的基于幾何方法的自由表面粒子判定，將該方法擴(kuò)展到了三維，并進(jìn)行了三維的球與平板碰撞模型計(jì)算，取得了良好的結(jié)果，但計(jì)算過程復(fù)雜，計(jì)算時(shí)間較長；Shibata等[10]對MPS方法開發(fā)新的壓力梯度模型時(shí)，通過光照法來判斷自由表面粒子，模擬容器內(nèi)靜止的流體時(shí)取得了良好的效果，但是這種方法不適用于運(yùn)動的流體。

Marrone等[11]對SPH方法中的自由表面粒子判定進(jìn)行了深入研究，首先通過求解重正化矩陣的特征值，初步篩選出自由表面粒子；然后通過幾何法，將第一步誤判的自由表面粒子剔除；使用該方法對潰壩模型進(jìn)行了模擬，取得了良好的結(jié)果，但是求解重正化矩陣較為復(fù)雜。也有一些學(xué)者通過修改計(jì)算參數(shù)來識別自由表面粒子。Lee等[4]在運(yùn)用MPS方法研究劇烈的自由表面振蕩問題時(shí)，通過對核函數(shù)、碰撞模型、源項(xiàng)以及梯度模型依次修正，達(dá)到了改進(jìn)判斷自由表面粒子的目的，并且減少了CPU的計(jì)算時(shí)間。本文針對MPS方法中自由表面粒子判定不準(zhǔn)確問題，開發(fā)了一種融合體積法和幾何法的自由表面粒子識別方法，并且對典型的潰壩問題進(jìn)行了數(shù)值模擬，進(jìn)而有效地改善了MPS方法中自由表面粒子的識別。

2 數(shù)值計(jì)算方法

2.1 控制方程

對于不可壓縮流體，其連續(xù)性方程和Navier-Stokes方程可表示為

(1)

(2)

2.2 傳統(tǒng)MPS方法

2.2.1 核函數(shù)

在MPS方法中，粒子間的相互作用關(guān)系通過核函數(shù)體現(xiàn)。在MPS方法中，核函數(shù)的公式如下，

(3)

式中r=|ri-rj|為兩個(gè)粒子之間的距離，ri為i粒子的位置，rj為i粒子支持域內(nèi)j粒子的位置，re為粒子作用域的半徑，一般re=2.1l0，其中l(wèi)0為初始粒子間距。

2.2.2 梯度模型

(4)

式中d為空間維度，n0為初始粒子數(shù)密度值，wi j=w(|rj-ri|)。在MPS法的模擬中，只有作為標(biāo)量的壓力值P參與上式計(jì)算。為了計(jì)算的穩(wěn)定性，MPS方法的壓力梯度計(jì)算一般為

(5)

Pj為i粒子支持域內(nèi)粒子的壓力，Pi min為i粒子支持域內(nèi)所有粒子中壓力最小的粒子壓力。

2.2.3 Laplace模型

(6)

式中d為空間維度，λ定義為

(7)

λ的引入是對有限范圍內(nèi)的核函數(shù)代替無限范圍的高斯函數(shù)帶來誤差的一種補(bǔ)償。Laplace模型用來離散二階導(dǎo)數(shù)項(xiàng)，在MPS法的模擬中，粘性力項(xiàng)用Laplace微分算子進(jìn)行離散。

2.2.4 壓力泊松方程

在傳統(tǒng)的MPS方法中，壓力通過求解壓力泊松方程得到，其中，系數(shù)矩陣項(xiàng)通過Laplace模型進(jìn)行離散，源項(xiàng)采用粒子數(shù)密度的偏差，可表示為

(8)

2.3 算法流程

MPS方法中，第一步計(jì)算重力和粘滯力作用下的速度增量，然后根據(jù)該速度增量對粒子的速度和位置進(jìn)行第一次的顯式修正;接下來計(jì)算此時(shí)的粒子數(shù)密度，建立并求解壓力泊松方程，然后計(jì)算出第二次的速度隱式修正；最后,根據(jù)此速度隱式修正粒子最終的速度和位置。MPS方法具體的模擬流程如圖1所示。

3 自由表面的判定

3.1 基于粒子數(shù)密度的自由表面識別方案

在最初的MPS方法中,流體的自由表面條件描述為，由于在自由表面以外沒有流體粒子，所以在自由表面處的流體粒子數(shù)密度小于流體內(nèi)部的粒子數(shù)密度，則采用式(9)來確定自由表面粒子，

(9)

3.2 改進(jìn)的自由表面識別方案

分步對自由表面粒子進(jìn)行篩選，首先采用幾何法對自由表面粒子進(jìn)行初步篩選；然后運(yùn)用體積法對自由表面粒子進(jìn)行進(jìn)一步識別，最終達(dá)到較好的自由表面粒子識別效果。

3.2.1 幾何法

(1) 基本思想

幾何法的基本思想是先以一個(gè)待識別的粒子為中心，畫一個(gè)直徑為h的圓；然后畫出以周圍粒子為中心的直徑為h的圓；如果這個(gè)待識別的自由表面粒子圓形未受到其周圍粒子的圓形全部覆蓋，則判定這個(gè)待識別的粒子為自由表面粒子。

圖2中灰色的粒子是需要判定的粒子，灰色圓的直徑為h；黑色的粒子是灰色粒子的周圍粒子，黑色圓形的半徑為h。假如以黑色粒子為中心的圓將以灰色粒子為中心的圓全部覆蓋，則灰色的粒子識別為內(nèi)部粒子；假如黑色圓形未將灰色圓形全部覆蓋，那么灰色的粒子判定為自由表面粒子。

(2) 實(shí)現(xiàn)方法

幾何法中，弧長的準(zhǔn)確計(jì)算是幾何法是否能夠準(zhǔn)確判定自由表面粒子的關(guān)鍵。如圖3所示，將待識別粒子的灰色圓心和周圍粒子的黑色圓心連接，那么兩點(diǎn)之間的距離為

(10)

式中 dx為兩點(diǎn)在x方向上的距離，dy為兩點(diǎn)在y方向上的距離。

圖1 MPS方法詳細(xì)模擬流程

Fig.1 Simulated flow chart of MPS method

∠B=arctan(dy/dx)

(11)

兩圓的圓心連線平分交點(diǎn)連線，如圖3所示，∠A是灰色圓心與兩圓交點(diǎn)的連接線和兩圓圓心的連線所成的夾角，

∠A=arccos(dr/h)

(12)

式中h為圓的直徑。圖3中，圓弧的起始角∠Start和終止角∠Stop的計(jì)算如下，

∠Start=∠B-∠A

(13)

∠Stop=∠B+∠A

(14)

實(shí)際上，灰色圓心與黑色圓心的連接線所處的象限不同，起始角和終點(diǎn)角的計(jì)算也不同。但是其基本方法和上述情況類似，限于文章篇幅不做具體展開。

在計(jì)算出黑色圓形截取灰色圓形的每段弧的起始角和終止角后，將這些弧依次排序，然后觀察其是否將灰色圓形全部覆蓋。如果全部覆蓋，則判定為內(nèi)部粒子；否則，判定為自由表面粒子,如圖5所示。

3.2.2 體積法

第一步已經(jīng)判斷出大部分自由表面粒子，但仍會出現(xiàn)個(gè)別內(nèi)部粒子誤判為自由表面粒子的情況。所以在第二步中，需要排除這些誤判的自由表面粒子。方法如下,假如一個(gè)自由表面粒子在其支持域內(nèi)的粒子都是流體粒子，則判定這個(gè)自由表面粒子為誤判的自由表面粒子。如圖5所示，黑色粒子識別為自由表面粒子，灰色粒子識別為流體粒子，黑色圓是黑色粒子的支持域?？梢钥闯?，在該黑色粒子的支持域內(nèi)的粒子都是內(nèi)部流體粒子，該黑色粒子為誤判粒子。由于該粒子在流體內(nèi)部，一定滿足

圖2 幾何法判定自由表面粒子示意圖

Fig.2 Surface particles detection of geometrical method

圖3 幾何法中弧長計(jì)算的示意圖

Fig.3 Arc calculation of geometrical method

(15)

4 數(shù)值模擬結(jié)果

4.1 計(jì)算模型及參數(shù)

圖4 幾何法中所截的不同弧排序的示意圖

Fig.4 Rank ordering of interceptive arc in geometrical method

圖5 幾何法中誤判為自由表面粒子的流體內(nèi)部粒子示意圖

Fig.5 Misjudgement of surface particles in geometrical method of the dam-break problem

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

圖7是三種自由表面粒子識別方案在t=0.4 s和t=1.1 s兩個(gè)時(shí)刻的自由表面粒子識別結(jié)果。此時(shí)流體流動較為平穩(wěn)，三種方案計(jì)算的潰壩流型是基本一致的，但是三種方法的自由表面判別結(jié)果卻有著明顯的差異。在t=0.4 s時(shí)刻，采用方案1進(jìn)行識別時(shí)，潰壩流體右上角處內(nèi)部有一層流體誤判為自由表面粒子，潰壩流體前鋒處有自由表面粒子未識別的情況。這是因?yàn)椴捎昧Ｗ訑?shù)密度大小進(jìn)行判定時(shí)，在模擬過程中，會造成表面粒子以及內(nèi)部粒子分布的不均勻，在粒子分布較為稀疏的地方，流體粒子的粒子數(shù)密度很可能小于判別系數(shù)β，那么此處的粒子就可能誤判為自由表面粒子。然而，采用幾何法自由表面粒子判定條件和采用幾何法+體積法自由表面粒子判定條件進(jìn)行自由表面粒子判定時(shí)，在潰壩流體右上角處的自由表面粒子都能夠準(zhǔn)確識別。在t=0.4 s時(shí)刻，采用方案2和方案3進(jìn)行識別時(shí)，潰壩流體右上角處的自由表面粒子都能準(zhǔn)確識別，在潰壩流體前鋒處也都沒有發(fā)生自由表面粒子未識別的情況。在t=1.1 s時(shí)刻，采用方案1進(jìn)行識別時(shí)，潰壩流體靠近前鋒處的流體出現(xiàn)了大量內(nèi)部流體誤判為自由表面粒子的情況，大概有20個(gè)流體內(nèi)部粒子誤判為自由表面粒子。采用方案2進(jìn)行識別時(shí)，潰壩流體靠近前鋒處的流體仍然有1個(gè)內(nèi)部流體粒子誤判為自由表面粒子。這是由于在流體流動過程中，粒子的流動極其復(fù)雜，在計(jì)算過程中仍然會出現(xiàn)幾何法未考慮到的特殊情況。針對這些特例，本文又增加了一次判斷，采用方案3進(jìn)行自由表面粒子識別時(shí)，就能避免這種特例的發(fā)生，潰壩流體靠近前鋒處再也沒有出現(xiàn)流體內(nèi)部粒子誤判為自由表面粒子的情況。

表1 不同自由表面識別方案

Tab.1 Different free surface detection methods

方案1 基于粒子數(shù)密度的識別方案方案2 基于幾何法的識別方案方案3 基于幾何法+體積法的識別方案

圖6 潰壩計(jì)算模型示意圖

Fig.6 Sketch of the dam-break problem

圖8是三種自由表面判定情況在t=2.45 s時(shí)刻的自由表面粒子分布示意圖，此時(shí)流體撞擊壁面形成了卷起的波浪，流體運(yùn)動非常劇烈。針對這種情況分析三種不同的自由表面粒子識別方案的效果。從圖9可以看出，當(dāng)采用方案1進(jìn)行判別時(shí)，由于卷起的波浪內(nèi)部流體運(yùn)動劇烈，導(dǎo)致粒子分布也非常不均勻，從而在粒子分布較稀疏處有很多流體內(nèi)部粒子誤判為自由表面粒子。所以方案1對流體運(yùn)動劇烈情況下的自由表面粒子識別效果不是很好。采用方案2進(jìn)行判別時(shí)，雖然卷起的波浪內(nèi)部流體粒子分布仍然不均勻，但是幾何法能夠較好地克服這一問題，之前誤判的大部分波浪內(nèi)部流體粒子都準(zhǔn)確地判定為自由表面粒子，說明方案2對自由表面粒子的識別效果優(yōu)于方案1。但可以看出，采用方案2時(shí)，卷起波浪內(nèi)部流體仍有個(gè)別粒子誤判為自由表面粒子。采用方案3進(jìn)行自由表面粒子判定時(shí)，所有流體內(nèi)部的粒子都沒有出現(xiàn)誤判為自由表面粒子的情況，對自由表面粒子的識別效果在三種方案中最優(yōu)。

4.3 計(jì)算效率的對比

圖9是三種自由表面判定方案在粒子數(shù)量分別為3772,4692,8348,11228,16762和24000下的計(jì)算時(shí)間，其中CPU 計(jì)算時(shí)間是指每循環(huán)10個(gè)時(shí)間步所需的時(shí)間?？梢钥闯觯S著粒子數(shù)量的增加，三種方案的CPU計(jì)算時(shí)間都會增加。其中，采用方案1所需的計(jì)算時(shí)間最少，采用方案2所需的計(jì)算時(shí)間次之，采用方案3所花費(fèi)的計(jì)算時(shí)間最多。這是因?yàn)榉桨?的自由表面識別最為簡單，只需要對比粒子數(shù)密度；而采用方案2時(shí)，判斷自由表面需要計(jì)算每個(gè)粒子受覆蓋的弧長，然后進(jìn)行排序，所以花費(fèi)的計(jì)算時(shí)間較多；采用方案3時(shí)，不但需要計(jì)算每個(gè)粒子受覆蓋的弧長且進(jìn)行排序，還需要對可能誤判的粒子進(jìn)行進(jìn)一步的排查，所以花費(fèi)的計(jì)算時(shí)間最多，但是自由表面的識別效果也最好?？梢钥闯?，隨著自由表面識別精度的提高，CPU計(jì)算時(shí)間也會相應(yīng)增加。雖然新的幾何法+體積法的自由表面識別方案在一定程度上降低了計(jì)算效率，但是降低的幅度不大。

圖7 三種自由表面判定情況在t=0.4 s和t=1.1 s時(shí)刻的自由表面粒子示意圖

Fig.7 Detection results of three different surface particles detection methods att=0.4 s andt=1.1 s

圖8 三種自由表面判定情況在t=2.15 s時(shí)刻的自由表面粒子示意圖

Fig.8 Detection results of three different surface particles detection methods att=2.45 s

圖9 三種自由表面識別方案在不同粒子數(shù)量下的模擬時(shí)間對比

Fig.9 Comparison of computational time between three different surface detection methods

5 結(jié) 論

本文采用MPS方法針對基于粒子數(shù)密度判定自由表面粒子效果不佳的問題，建立了兩步法，包括幾何法和體積法進(jìn)行自由表面粒子判定的方法。采用新的兩步法自由表面粒子判定方法，對潰壩模型問題進(jìn)行了數(shù)值模擬。并且對潰壩流體在不同的流體運(yùn)動劇烈程度下，采用不同的自由表面粒子判定方法時(shí)，比較分析自由表面粒子判定效果。結(jié)果表明，本文新建立的兩步法判定自由表面粒子對于潰壩流體流動時(shí)能夠準(zhǔn)確地判斷其自由表面粒子，克服了之前基于粒子數(shù)密度判定自由表面粒子的缺陷，為以后采用MPS方法研究兩相流問題，不同流體之間通過界面進(jìn)行傳熱傳質(zhì)問題打下了良好的基礎(chǔ)。