， ， ，

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

1 引 言

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

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

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

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

2 數(shù)值計算方法

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|為兩個粒子之間的距離，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參與上式計算。為了計算的穩(wěn)定性，MPS方法的壓力梯度計算一般為

(5)

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

2.2.3 Laplace模型

(6)

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

(7)

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

2.2.4 壓力泊松方程

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

(8)

2.3 算法流程

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

3 自由表面的判定

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

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

(9)

3.2 改進的自由表面識別方案

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

3.2.1 幾何法

(1) 基本思想

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

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

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

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

(10)

式中 dx為兩點在x方向上的距離，dy為兩點在y方向上的距離。

圖1 MPS方法詳細模擬流程

Fig.1 Simulated flow chart of MPS method

∠B=arctan(dy/dx)

(11)

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

∠A=arccos(dr/h)

(12)

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

∠Start=∠B-∠A

(13)

∠Stop=∠B+∠A

(14)

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

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

3.2.2 體積法

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

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

Fig.2 Surface particles detection of geometrical method

圖3 幾何法中弧長計算的示意圖

Fig.3 Arc calculation of geometrical method

(15)

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

4.1 計算模型及參數(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 計算結(jié)果及分析

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

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

Tab.1 Different free surface detection methods

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

圖6 潰壩計算模型示意圖

Fig.6 Sketch of the dam-break problem

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

4.3 計算效率的對比

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

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

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

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

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

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

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

5 結(jié) 論

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