采用多邊形壁面的MPS粒子分裂算法研究

張 靜，熊進(jìn)標(biāo)

(上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240)

反應(yīng)堆嚴(yán)重事故模擬中涉及大量含有自由界面或組分界面的多相流動(dòng)，不具有拓?fù)浣Y(jié)構(gòu)的粒子法在這類流動(dòng)的模擬中有其獨(dú)特的優(yōu)勢(shì)。李勇霖等[1]基于移動(dòng)粒子半隱式(MPS)法[2]模擬堆芯材料在高溫熔化過程中的共晶反應(yīng)現(xiàn)象。張明昊等[3]采用MPS-MAFL方法針對(duì)入口流量脈動(dòng)條件下的單氣泡垂直上升運(yùn)動(dòng)行為進(jìn)行了特性分析。Zhu等[4-6]引入固體被動(dòng)漂移模型、勢(shì)能界面張力等模型驗(yàn)證MPS法對(duì)熔融現(xiàn)象計(jì)算的可靠性。但由于傳統(tǒng)的MPS法中粒子尺寸相等，無法對(duì)局部區(qū)域(如相變界面附近)進(jìn)行加密計(jì)算。因此提出多粒徑模型，Tang等[7]針對(duì)不同分辨率界面上的粒子，對(duì)梯度算子模型和拉普拉斯算子模型添加修正系數(shù)來滿足守恒性。Chen等[8]采用三次樣條函數(shù)代替原先的核函數(shù)，添加額外的權(quán)重以及采取5步分裂方法來保證分裂計(jì)算的準(zhǔn)確性和壓力的穩(wěn)定性。Shibata等[9]設(shè)置局部重疊的子區(qū)域，并添加區(qū)域邊界條件來實(shí)現(xiàn)局部的精細(xì)化計(jì)算。

此外，為了改善壁面粒子對(duì)計(jì)算效率和準(zhǔn)確性產(chǎn)生的影響，人們提出了多邊形壁面模型。與傳統(tǒng)的采用固定流體粒子作為壁面的方式相比，該模型不設(shè)置壁面粒子，因此在求解壓力泊松方程時(shí)僅包含流體粒子，提高迭代計(jì)算效率；同時(shí)，壁面的壓力采取合理的方式替代求解。對(duì)于壁面附近粒子數(shù)密度以及壓力的計(jì)算方法，Harada等[10]基于粒子在壁面壓力梯度的作用下回到平衡位置的假設(shè)提出多邊形壁面形式。Zhang等[11]在Harada的方法上進(jìn)行改進(jìn)，提出一種新的虛擬壁面粒子生成方式，改進(jìn)適用于多邊形壁面的壓力源項(xiàng)[12]。Mitsume等[13]通過計(jì)算鏡面粒子的壓力梯度代替實(shí)際壁面梯度。

在目前的研究進(jìn)展中，粒子分裂一般需要設(shè)置較為復(fù)雜的邊界條件，采取多邊形壁面模型的算法中往往設(shè)置壓力顯示計(jì)算的方式，鮮有粒子分裂模型與多邊形壁面模型相結(jié)合的研究。本文針對(duì)MPS進(jìn)行算法改進(jìn)，使其能適用于多邊形壁面條件下基于壓力隱式求解的粒子分裂計(jì)算，對(duì)比潰壩實(shí)驗(yàn)驗(yàn)證改進(jìn)后模型的準(zhǔn)確性和高效性，為進(jìn)一步計(jì)算三維相變傳熱模型奠定基礎(chǔ)。

1 移動(dòng)粒子半隱式法

對(duì)于不涉及傳熱的流動(dòng)過程，MPS的控制方程包括連續(xù)性方程和動(dòng)量方程：

(1)

(2)

式中：ρ為密度；t為時(shí)間；u為速度；P為壓力；ν為運(yùn)動(dòng)黏性；g為重力加速度；f為重力以外的其他外力。

MPS法將求解域離散為一系列的粒子，通過粒子間相互作用模型離散控制方程，基于追蹤粒子得到計(jì)算域的信息。粒子相互作用模型包括核函數(shù)、梯度算子模型和拉普拉斯算子模型3個(gè)基本模型，常用的核函數(shù)模型為：

re=RDi

(3)

式中：|rij|為粒子i與j的間距；re為有效半徑；R為有效半徑系數(shù)，通常取2.1～4.1；Di為粒子直徑。

對(duì)有效半徑內(nèi)粒子的核函數(shù)求和，得到目標(biāo)粒子的粒子數(shù)密度ni為：

ni=∑wi(|rij|)

(4)

梯度算子模型為：

(5)

拉普拉斯算子模型為：

(6)

通過聯(lián)立動(dòng)量和連續(xù)性方程，可得壓力泊松方程為：

(7)

在壓力計(jì)算中，表面粒子的壓力設(shè)置為0。通常采用粒子數(shù)密度閾值來判定表面粒子[14]，即：

〈n〉i<αn0

(8)

MPS算法流程如圖1所示。在1個(gè)時(shí)間步長內(nèi)，首先利用黏性力與重力顯式更新速度和位置，然后通過式(7)隱式計(jì)算壓力，并利用式(5)得到壓力梯度，再次更新粒子的速度及位置信息，確認(rèn)達(dá)到收斂標(biāo)準(zhǔn)后進(jìn)入下一時(shí)間步的計(jì)算。

圖1 MPS算法流程

2 模型開發(fā)

在多粒徑計(jì)算中，直徑不同的粒子質(zhì)量存在差異，傳統(tǒng)MPS法的粒子數(shù)密度計(jì)算方法無法準(zhǔn)確反映粒子聚集程度(即密度變化)，從而導(dǎo)致壓力計(jì)算不準(zhǔn)確。同時(shí)大量的壁面粒子也增大了泊松方程的求解量。本文針對(duì)上述問題進(jìn)行改進(jìn)，開發(fā)適用于多邊形壁面的粒子分裂模型。

2.1 有效半徑

當(dāng)計(jì)算域內(nèi)的粒子尺寸不等時(shí)，采用傳統(tǒng)的有效半徑計(jì)算方法會(huì)影響計(jì)算的準(zhǔn)確性。本研究中僅考慮兩種粒徑的粒子計(jì)算，當(dāng)無外力作用下粒子均勻排布時(shí)，如圖2a所示，粒徑不同的i粒子與j粒子位于分辨率界面上，實(shí)線圈為i粒子的作用域，虛線圈為j粒子的作用域。理論上兩者的粒子數(shù)密度應(yīng)相等，粒子間保持相對(duì)靜止；而由式(3)、(4)計(jì)算得到的i粒子的數(shù)密度遠(yuǎn)大于j粒子的數(shù)密度，兩個(gè)粒子間會(huì)生成較大的作用力而運(yùn)動(dòng)，因此需要對(duì)有效半徑進(jìn)行改進(jìn)。

本文參考Tanaka等[15]的處理方法，采用平均有效半徑代替原固定的有效半徑，即：

(9)

式中，R取值為3.1。

改進(jìn)后的作用域如圖2b所示。當(dāng)大粒子i與小粒子相互作用時(shí)，作用域如小半圓實(shí)線圈所示，作用范圍與改進(jìn)前相比有所縮小；粒子i與其直徑相等的大粒子作用時(shí)，作用域如大半圓實(shí)線圈所示，與原有效半徑相等。同理，小粒子j的作用域也改進(jìn)為大小半圓組合的不規(guī)則范圍。改進(jìn)后的有效半徑計(jì)算方法減小了粒子i與粒子j間的數(shù)密度差值，同時(shí)也避免了由于兩個(gè)粒子的有效半徑的不同而可能造成的相互作用力不相等現(xiàn)象。

a——改進(jìn)前；b——改進(jìn)后

2.2 粒子數(shù)密度

粒子數(shù)密度ni包括流體粒子的密度nif及固體壁面的密度niw：

ni=nif+niw

(10)

nif的計(jì)算參考多粒徑模型的處理方法[7]，根據(jù)牛頓第三定律對(duì)數(shù)密度添加系數(shù)項(xiàng)：

(11)

式中，m為粒子質(zhì)量。

niw的求解與傳統(tǒng)MPS法計(jì)算不同。首先在初始化過程中，對(duì)計(jì)算域建立大小均勻的背景網(wǎng)格，網(wǎng)格尺寸與小粒子尺寸一致。計(jì)算每個(gè)網(wǎng)格點(diǎn)到多邊形壁面的最小距離，并在多邊形壁面外側(cè)生成虛擬壁面粒子。對(duì)位于壁面有效半徑內(nèi)的網(wǎng)格點(diǎn)，計(jì)算其曲率系數(shù)Cg[12]為：

(12)

式中：Zg為虛擬壁面粒子在g網(wǎng)格點(diǎn)處的數(shù)密度；Wg為同等距離下，固體粒子組成平板壁面時(shí)對(duì)應(yīng)的g網(wǎng)格處的壁面數(shù)密度。當(dāng)Cg=1時(shí)，在g網(wǎng)格點(diǎn)有效半徑內(nèi)壁面為平板；當(dāng)Cg≠1時(shí)，該網(wǎng)格點(diǎn)有效半徑內(nèi)存在壁面轉(zhuǎn)角。

在二維計(jì)算中，搜索流體粒子周圍的4個(gè)網(wǎng)格點(diǎn)，并基于網(wǎng)格點(diǎn)的壁面距離插值得到流體粒子的壁面距離riw，進(jìn)而計(jì)算該距離對(duì)應(yīng)的平板壁面數(shù)密度Wiw；對(duì)網(wǎng)格的曲率系數(shù)進(jìn)行插值得到粒子在該位置處的曲率系數(shù)Ciw，將平板壁面數(shù)密度Wiw與曲率系數(shù)Ciw相乘，得到流體粒子的實(shí)際壁面數(shù)密度niw。壁面數(shù)密度計(jì)算示意圖如圖3所示，紅點(diǎn)為參與計(jì)算的網(wǎng)格點(diǎn)，綠線為多邊形壁面，虛線粒子為虛擬壁面粒子。

圖3 壁面數(shù)密度計(jì)算示意圖

2.3 黏性拉普拉斯算子模型

(13)

(14)

(15)

式中，uw為壁面速度。

2.4 壓力泊松方程

采用多邊形壁面模型后，無需設(shè)置壁面粒子。因此在求解壓力泊松方程時(shí)，僅考慮流體粒子間的相互作用。多粒徑的壓力拉普拉斯算子模型的計(jì)算采取式(16)代替式(6)：

(16)

參考Zhang等[12]對(duì)適用于多邊形壁面的壓力源項(xiàng)的處理，添加粒子流體數(shù)密度與粒子總數(shù)密度的比值這一系數(shù)，來反映壁面壓力對(duì)流體的影響。本文將其改進(jìn)為適用于分裂計(jì)算的表達(dá)形式：

(17)

式中：γ為修正系數(shù)，取0.015；k為當(dāng)前時(shí)間步，k+1為下一時(shí)間步；u*為兩個(gè)時(shí)間步間的臨時(shí)速度。聯(lián)立式(16)與(17)，隱式求解得到流體粒子的壓力。

2.5 梯度算子模型

(18)

(19)

圖4 鏡面粒子壓力梯度計(jì)算[16]

(20)

2.6 自由表面粒子識(shí)別

為了提高自由表面粒子識(shí)別的準(zhǔn)確性，本文參考Shibata等[17]提出的基于弧度的表面識(shí)別方法進(jìn)行模型改進(jìn)。模型中，將目標(biāo)粒子i假想為一個(gè)點(diǎn)光源，照亮半徑為2.1di的圓周，計(jì)算有效半徑內(nèi)的相鄰粒子由于遮擋光線而在圓周上形成的陰影弧長。定義表面率A為：

(21)

由于采用式(21)進(jìn)行表面識(shí)別需要對(duì)所有粒子進(jìn)行判斷，計(jì)算量較大，本文采用式(8)與式(21)相結(jié)合的方式。首先對(duì)所有粒子進(jìn)行基于粒子數(shù)密度的篩選，α取0.97，將符合條件的粒子再進(jìn)行基于弧度的識(shí)別判斷，從而在保證精度的同時(shí)提高效率。

2.7 分裂模型

粒子分裂計(jì)算的流程如圖5所示。

圖5 粒子分裂計(jì)算流程圖

首先設(shè)置分裂條件，選取粒子位置為判斷標(biāo)準(zhǔn)，在計(jì)算域內(nèi)劃分出需要精細(xì)計(jì)算的分裂區(qū)域。識(shí)別出符合分裂條件的大粒子，對(duì)分裂生成新的小粒子進(jìn)行位置賦值。本文僅考慮1個(gè)大粒子均勻分裂成4個(gè)小粒子的情況，其排列分布如圖6所示。為了保證分裂后粒子速度與周圍粒子光滑過渡，本文參考Tanaka等[15]提出的梯度修正法，給定新粒子的速度：

圖6 粒子分裂后位置及大小示意圖

(22)

式中：ui,M為分裂后粒子速度；ui為分裂前速度；ri,M為分裂后粒子位置；ri為分裂前位置；M為分裂后粒子序號(hào)，M=1、2、3、4。

在對(duì)體積等常量進(jìn)行賦值后，將發(fā)生分裂的大粒子的類型更新為Ghost，不參與后續(xù)的計(jì)算。分裂循環(huán)完成后，更新全計(jì)算域粒子的數(shù)密度，結(jié)束分裂程序的計(jì)算。

3 潰壩模擬

3.1 靜壓測(cè)試

靜壓計(jì)算中理論上粒子處于靜止?fàn)顟B(tài)，監(jiān)測(cè)點(diǎn)的壓力具有精確的理論值，因此采用靜壓計(jì)算來評(píng)價(jià)改進(jìn)后MPS模型的準(zhǔn)確性以及穩(wěn)定性。液體幾何長度為0.36 m、高度為0.48 m，坐標(biāo)原點(diǎn)位于左下角，監(jiān)測(cè)(0.2 m，0.01 m)點(diǎn)處的壓力。計(jì)算多邊形壁面條件下，粒子直徑分別為0.01 m及0.005 m工況(背景網(wǎng)格尺寸均為0.005 m)下的壓力波動(dòng)，如圖7所示。計(jì)算達(dá)到穩(wěn)定后，兩個(gè)工況下壓力均波動(dòng)較小且略低于理論值，相對(duì)誤差約為6.25%。因此，靜壓測(cè)試可驗(yàn)證多邊形壁面模型的準(zhǔn)確性和穩(wěn)定性。

圖7 不同工況在監(jiān)測(cè)點(diǎn)的壓力

3.2 無擋板潰壩

計(jì)算Martin等[18]的無擋板潰壩實(shí)驗(yàn)，驗(yàn)證粒子分裂模型的可靠性。無擋板潰壩模型如圖8所示，水柱高0.44 m、長0.22 m。設(shè)置粒子的密度為1 000 kg/m3，運(yùn)動(dòng)黏性為1.0 mm2/s，重力加速度為9.81 m/s2，時(shí)間步長為1×10-4s，共計(jì)算了3種工況：1)粒子直徑為0.01 m；2)粒子直徑為0.005 m；3)粒子在波前距離大于0.4 m時(shí)分裂，分裂前直徑為0.01 m、分裂后直徑為0.005 m，分裂后粒子排布形式如圖6所示。

圖8 無擋板潰壩模型

圖9 無量綱時(shí)間變化圖

3.3 有擋板潰壩

計(jì)算有檔板的潰壩模型[19]，以驗(yàn)證壓力計(jì)算的準(zhǔn)確性以及分裂模型的計(jì)算效率。有擋板潰壩模型如圖10所示。水柱高0.12 m、長0.68 m，容器長1.18 m，壓力測(cè)試點(diǎn)位于右側(cè)壁面距離地面0.01 m處，即A點(diǎn)所示。

圖10 有擋板潰壩模型

計(jì)算了3種工況：1)壁面設(shè)置為傳統(tǒng)固定粒子方式，粒子直徑為0.01 m；2)壁面設(shè)置為傳統(tǒng)固定粒子方式，粒子直徑為0.005 m；3)壁面設(shè)置為多邊形壁面形式，當(dāng)粒子的橫坐標(biāo)(長度方向)大于0.7 m時(shí)分裂，分裂前粒子直徑為0.01 m，分裂后粒子直徑為0.005 m，分裂后粒子排布形式如圖6所示。設(shè)置時(shí)間步長為5×10-4s，選取0.5～1.0 s共6個(gè)時(shí)刻的壓力計(jì)算結(jié)果進(jìn)行比較，如圖11所示?？煽闯觯煌瑫r(shí)刻的3種工況下的粒子流動(dòng)具有相似性。對(duì)于自由表面的捕捉，采用均勻小粒子計(jì)算以及分裂計(jì)算的工況(圖11b、c)能得到清晰的液面，而均勻大粒子的工況(圖11a)得到的結(jié)果較為粗糙。

a——粒子壁面，均勻粒徑0.01 m；b——粒子壁面，均勻粒徑0.005 m；c——多邊形壁面，初始粒徑0.01 m，分裂后0.005 m

在計(jì)算過程中對(duì)圖10中的A點(diǎn)進(jìn)行壓力檢測(cè)，并與實(shí)驗(yàn)值進(jìn)行比較(圖12)。無分裂+粒子壁面工況的結(jié)果為藍(lán)色線所示，撞擊壁面時(shí)間、撞擊后監(jiān)測(cè)點(diǎn)的壓力均與實(shí)驗(yàn)值吻合良好；分裂+粒子壁面工況的結(jié)果如紅色線所示，除壓力次高峰與實(shí)驗(yàn)值相比滯后約0.1 s外，其余均與實(shí)驗(yàn)值吻合良好；分裂+多邊形壁面工況的結(jié)果為綠線所示，撞擊壁面時(shí)間與實(shí)驗(yàn)值相近，但是檢測(cè)點(diǎn)的壓力波動(dòng)幅值較大，且壓力偏大。因此改進(jìn)后模型的壓力穩(wěn)定性需要進(jìn)一步的研究和改善。

圖12 A點(diǎn)壓力比較

對(duì)改進(jìn)后模型的計(jì)算效率進(jìn)行分析，計(jì)算均采用i9-10900處理器，內(nèi)核頻率為2.8 GHz。對(duì)數(shù)密度計(jì)算模塊、壓力求解模塊以及計(jì)算總時(shí)長進(jìn)行比較(表1)。其中，粒子個(gè)數(shù)為整個(gè)計(jì)算過程中的平均粒子個(gè)數(shù)，模塊計(jì)算時(shí)長為1個(gè)時(shí)間步內(nèi)、單次計(jì)算所需要的平均時(shí)間；總時(shí)長占比為改進(jìn)后的模型與傳統(tǒng)MPS模型(模型1)總時(shí)長的比值。由表1可看出，若僅采用分裂計(jì)算(模型1與2比較)，粒徑增大、粒子個(gè)數(shù)減少使得數(shù)密度及壓力求解模塊的時(shí)長減小，總計(jì)算時(shí)長為原來的44.28%；若僅采用多邊形壁面(模型1與3比較)，參與壓力泊松求解的粒子數(shù)減少、迭代效率提高，總時(shí)長占比為37.42%；若采用多邊形壁面與分裂模型相結(jié)合的形式(模型1與4比較)，所需粒子個(gè)數(shù)進(jìn)一步減少，平均僅需1 486.3個(gè)，總時(shí)長占比為18.75%，略高于兩個(gè)模型單獨(dú)作用的時(shí)長占比的乘積(16.57%)。因此，采用多邊形壁面的粒子分裂模型可大幅提高計(jì)算的效率。

表1 不同工況計(jì)算時(shí)長比較

4 結(jié)論

本文基于MPS法，開發(fā)適用于多邊形壁面的粒子分裂模型，將改進(jìn)后的模型進(jìn)行無擋板以及含有擋板的潰壩實(shí)驗(yàn)計(jì)算，得出以下結(jié)論：1)改進(jìn)后的模型流動(dòng)結(jié)果與實(shí)驗(yàn)值吻合良好；2)改進(jìn)后的模型壓力波動(dòng)較大且數(shù)值偏高，需要對(duì)壓力穩(wěn)定性進(jìn)行改善；3)采用多邊形壁面的粒子分裂模型可大幅提高M(jìn)PS法的計(jì)算效率，改進(jìn)后模型計(jì)算總時(shí)長與原模型的占比為18.75%。