張 靜，熊進標

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

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

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

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

1 移動粒子半隱式法

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

(1)

(2)

式中：ρ為密度；t為時間；u為速度；P為壓力；ν為運動黏性；g為重力加速度；f為重力以外的其他外力。

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

re=RDi

(3)

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

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

ni=∑wi(|rij|)

(4)

梯度算子模型為：

(5)

拉普拉斯算子模型為：

(6)

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

(7)

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

〈n〉i<αn0

(8)

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

圖1 MPS算法流程

2 模型開發(fā)

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

2.1 有效半徑

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

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

(9)

式中，R取值為3.1。

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

a——改進前；b——改進后

2.2 粒子數(shù)密度

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

ni=nif+niw

(10)

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

(11)

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

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

(12)

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

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

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

2.3 黏性拉普拉斯算子模型

(13)

(14)

(15)

式中，uw為壁面速度。

2.4 壓力泊松方程

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

(16)

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

(17)

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

2.5 梯度算子模型

(18)

(19)

圖4 鏡面粒子壓力梯度計算[16]

(20)

2.6 自由表面粒子識別

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

(21)

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

2.7 分裂模型

粒子分裂計算的流程如圖5所示。

圖5 粒子分裂計算流程圖

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

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

(22)

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

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

3 潰壩模擬

3.1 靜壓測試

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

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

3.2 無擋板潰壩

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

圖8 無擋板潰壩模型

圖9 無量綱時間變化圖

3.3 有擋板潰壩

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

圖10 有擋板潰壩模型

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

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

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

圖12 A點壓力比較

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

表1 不同工況計算時長比較

4 結(jié)論

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