肖鑫坤，蔡慶航，陳榮華,*，丁 雯，張 魁，郭凱倫，田文喜，秋穗正，蘇光輝

(1.西安交通大學(xué) 動力工程多相流國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

移動粒子半隱式(MPS)方法[1]已有近30年的歷史，其主要思想是將流體用一系列離散的帶有拉格朗日屬性的粒子表示，將流體運(yùn)動控制方程的各項(xiàng)微分算子通過粒子間相互作用來表達(dá)，以實(shí)現(xiàn)對流體運(yùn)動的準(zhǔn)確模擬。相較于網(wǎng)格方法，MPS方法不存在網(wǎng)格畸變這一固有缺陷，而其拉格朗日特性使得MPS方法在捕捉自由表面和界面運(yùn)動、捕捉相態(tài)變化、固體變形等方面具有獨(dú)特的優(yōu)勢，因此得到的了眾多研究者的關(guān)注，并在核工程領(lǐng)域得到了廣泛應(yīng)用。

在核反應(yīng)堆嚴(yán)重事故分析領(lǐng)域，通過在MPS方法中引入傳熱相變模型、黏度變化模型[2]和共晶反應(yīng)模型，使用MPS方法對管內(nèi)流動凝固[3-5]、沸水堆燃料支撐件內(nèi)的熔融物流動凝固[6-7]、鉛鉍平板消熔[7]、共晶反應(yīng)[8-9]、熔融物與混凝土相互作用[10-11]、燃料元件熔化[12-13]、鋯水氧化反應(yīng)[14]和ADS顆粒靶熔化[15]等核反應(yīng)堆嚴(yán)重事故關(guān)鍵現(xiàn)象展開了模擬研究，這些研究驗(yàn)證了MPS方法在傳熱相變和物質(zhì)變化中應(yīng)用的可行性和準(zhǔn)確性，大大拓寬了MPS方法的應(yīng)用范圍。

在核反應(yīng)堆中發(fā)生堆芯熔化嚴(yán)重事故時，堆芯熔化后產(chǎn)生的熔融物發(fā)生遷徙，未發(fā)生熔化的燃料微粒、基體碎片、包殼碎片也被熔融物夾帶著發(fā)生遷徙，這是一個包含大量離散固體的流固耦合問題。已有的嚴(yán)重事故程序?qū)Χ研救刍瘑栴}進(jìn)行研究時，大都采用集總參數(shù)的燭化模型，如MELCOR等，這種模型計(jì)算得到的結(jié)果較為粗糙，且不考慮未熔化的離散固體之間的相互作用。

本研究使用MPS方法對這種現(xiàn)象進(jìn)行模擬，對熔融物夾帶的固體之間的相互作用以及液相與固相之間的作用進(jìn)行分析?，F(xiàn)有的處理方法難以進(jìn)行準(zhǔn)確的模擬，需要進(jìn)行相關(guān)模型的開發(fā)。

本文在基于MPS方法的核反應(yīng)堆關(guān)鍵熱工安全現(xiàn)象分析軟件平臺(PANDA)上進(jìn)一步發(fā)展離散單元法(DEM)耦合模塊(PANDA-DEM)，在使用的流固耦合模型中引入多相黏性的處理方法提高穩(wěn)定性和準(zhǔn)確性，并通過滑坡算例對本模塊使用的流固耦合模型的穩(wěn)定性和準(zhǔn)確性進(jìn)行驗(yàn)證。

1 數(shù)值模型

使用MPS方法對液相進(jìn)行描述，分析流體的運(yùn)動；使用DEM對固相進(jìn)行描述，計(jì)算固體粒子之間的相互作用力；液相和固相之間的相互作用力被歸納為黏性力和壓力作用。

1.1 液相

使用MPS方法來對流體的運(yùn)動進(jìn)行模擬，控制方程包括連續(xù)性方程和動量守恒方程：

(1)

(2)

式中：ρ為密度；t為時間；u為速度向量；p為壓力；μ為黏性系數(shù)；g為重力向量；f為表面張力。

流體的運(yùn)動過程通過上述控制方程進(jìn)行約束并得到，具體的思路是先顯式計(jì)算動量守恒方程中的黏性項(xiàng)、重力項(xiàng)和表面張力項(xiàng)，得到估算的速度和位置，通過不可壓縮流體的連續(xù)性方程建立壓力泊松方程，使用求解得到的壓力梯度對估算得到的粒子速度和位置進(jìn)行修正。

1.2 固相

使用DEM[16]求解離散固體的整體運(yùn)動狀態(tài)。在DEM理論中，離散的結(jié)構(gòu)單元被稱為顆粒，具有與MPS方法中的粒子相似的性質(zhì)。這兩種方法都是將介質(zhì)進(jìn)行離散化處理，通過對離散單元的運(yùn)動分析得到介質(zhì)整體的運(yùn)動狀態(tài)，因此能夠在同樣的框架內(nèi)進(jìn)行計(jì)算。

DEM的主要思路是對每一個離散單元進(jìn)行受力分析，通過牛頓第二定律得到其運(yùn)動狀態(tài)。

首先分析離散單元的運(yùn)動方程，對于當(dāng)前時刻，任意離散單元均可由牛頓第二定律得到兩個加速度：

(3)

(4)

下一時刻離散單元的速度u和角速度ω、位移r和旋轉(zhuǎn)角度θ為：

(5)

(6)

(7)

(8)

式中，ΔtDEM為DEM計(jì)算的時間步長。

由此可得到單個離散單元隨時間的運(yùn)動。離散單元加速度的合力與合力矩通過下式計(jì)算：

(9)

(10)

ΔTXe,j→i=krΔθX,i,j

(11)

Td,j→i=crΔθX,i,j/ΔtDEM

(12)

式中：e和d分別為彈性力和阻尼力；k和c分別為彈性系數(shù)和阻尼系數(shù)；下標(biāo)n和s分別代表接觸面的法向和切向，其中y軸方向和z軸方向同為切線方向，計(jì)算時只有方向不同，其余參數(shù)均相同；FX，j→i和FY(Z)，j→i為離散單元j對離散單元i的法向作用力和切向作用力；ΔXi，j和ΔY(Z)i，j為離散單元j和離散單元i在坐標(biāo)系x、y、z3軸方向上的相對位移；ΔTXe，j→i為x軸方向的回轉(zhuǎn)力矩的彈性增量；ΔθX，i，j為離散單元j與離散單元i的相對轉(zhuǎn)角；Td，j→i為離散單元j與離散單元i之間回轉(zhuǎn)力矩的阻尼分量；下標(biāo)r表示回轉(zhuǎn)。需要指出的是，計(jì)算固體之間的相互作用力時需要使用的相對位移通過速度和時間步長得到，而兩個固體粒子之間的相對位置通過上一個時間步長的位置得到。

1.3 流固耦合

通過將DEM計(jì)算離散固體之間的相互作用力視為MPS方法的一個模型來進(jìn)行二者之間的耦合，也就是在PANDA中開發(fā)DEM耦合模塊PANDA-DEM。將DEM的離散單元視為MPS方法中的粒子，二者使用相同的粒子直徑和作用半徑。對于離散固體體系受到的外力，也就是流體對固體的作用力，通過式(13)[17]進(jìn)行計(jì)算：

(13)

式中：下標(biāo)ls表示流固耦合界面；f為界面上的力。

式(13)認(rèn)為，流體對固體的作用力主要是流體對固體的壓力作用和流體對固體的黏性力作用，當(dāng)然需要認(rèn)識到這樣的作用只發(fā)生在流體與固體接觸的部分，也就是固體粒子之間的相互作用只需要進(jìn)行DEM計(jì)算，只對與流體粒子發(fā)生接觸的那部分固體粒子進(jìn)行MPS方法中的黏性項(xiàng)和壓力項(xiàng)計(jì)算。

此時，固體粒子對流體粒子的作用也就是流體粒子與固體粒子之間的耦合作用力fls的反向作用力，會被包含在流體粒子受到的黏性力和壓力中，因?yàn)榱黧w粒子的這兩個力是通過作用域內(nèi)全部粒子對它的作用力的積分得到，包括固體粒子。對于固體粒子而言，需要針對不同的情況進(jìn)行分析。對于存在于流固相界面上的固體粒子，其控制方程為：

(14)

對于不在相界面上的固體粒子：

(15)

式中：fDEM為使用DEM方法計(jì)算得到的固體粒子之間的相互作用力；ρs和ρl為固體的密度和流體的密度；ρsg為固體粒子受到的重力，對于fDEM而言，它是通過DEM模型計(jì)算得到的由于固體之間的相互作用產(chǎn)生的力，這個力通過兩個時層之間粒子的相對位移來計(jì)算得到。

在計(jì)算流體粒子與固體粒子之間的耦合作用力中的黏性項(xiàng)時，為充分考慮固體粒子與流體粒子之間存在的物性差異，體現(xiàn)固體粒子的物理性質(zhì)，也為了避免在流體與固體之間產(chǎn)生一個尖銳的相界面，通過引入多相流的處理思路，給予離散的固體一個與其物理性質(zhì)相關(guān)的黏性系數(shù)，并通過調(diào)和黏性系數(shù)來計(jì)算黏性項(xiàng)：

(16)

式中，μij為調(diào)和黏性系數(shù)，不難看出μij在μi和μj相等時退回單相，而兩個黏性系數(shù)不同只在相界面上出現(xiàn)，這是一個廣義模型。

此時黏性項(xiàng)采用式(17)進(jìn)行計(jì)算：

G(‖rj-ri‖,re)

(17)

通過計(jì)算發(fā)生接觸的流體粒子與固體粒子之間的黏性力與壓力作用，將這兩個力作為DEM計(jì)算時的外部輸入，就能得到離散的固體在流體作用下的運(yùn)動，也就實(shí)現(xiàn)了對流體與固體之間相互作用的耦合分析。

一般情況下，MPS方法計(jì)算的時間步長要大于DEM計(jì)算，時間步長的耦合通過先對流體粒子推進(jìn)一個MPS方法的時間步長，如果存在DEM計(jì)算，就進(jìn)入DEM計(jì)算內(nèi)循環(huán)，推進(jìn)數(shù)個DEM時間步長，令其總的時間步長等于MPS方法計(jì)算時間步長來實(shí)現(xiàn)。

PANDA-DEM的算法流程如下：1)建立研究對象的初始狀態(tài)模型，輸入初始參數(shù)；2)初始化速度和位置參數(shù)，用于后續(xù)的計(jì)算；3)推進(jìn)1個時間步長，使用MPS方法計(jì)算粒子的速度和位置；4)判斷是否發(fā)生固體粒子之間的相互碰撞，如果發(fā)生，進(jìn)入DEM計(jì)算內(nèi)循環(huán)，如果不發(fā)生，結(jié)束當(dāng)前時間步長，進(jìn)入下一時間步長；5)進(jìn)入DEM內(nèi)循環(huán)后，使用DEM模型計(jì)算固體粒子的速度和位置；6)結(jié)束內(nèi)循環(huán)，進(jìn)入下一時間步長。

2 滑坡算例的模擬

通過兩個滑坡算例對模型的穩(wěn)定性與準(zhǔn)確性進(jìn)行分析：算例1是固體滑坡算例，即純離散固體顆粒沿滑坡的下落過程，能夠驗(yàn)證耦合分析模型對于固體運(yùn)動的計(jì)算能力；算例2是水下滑坡算例，即離散的固體顆粒沿斜坡下落的過程中與流體發(fā)生耦合，流體流場與固體運(yùn)動都會受到流固耦合作用的影響。

2.1 固體滑坡算例

固體滑坡算例是干顆粒材料(即玻璃微珠)沿45°光滑斜面的滑動，圖1示出計(jì)算域的初始粒子排列。在實(shí)驗(yàn)開始時，將擋板快速移出，顆粒狀材料受到重力的作用沿滑坡向下滑動。固體滑坡算例基于Tajnesaie等[19]的實(shí)驗(yàn)設(shè)計(jì)得到，需要使用的幾何參數(shù)為：h1=0.2 m，h2=0.14 m，α=45°。DEM計(jì)算需要使用的參數(shù)為：法向彈性系數(shù)kE=2×105N/m，法向阻尼系數(shù)cE=30 N/m，切向彈性系數(shù)kN=2×104N/m，切向阻尼系數(shù)cN=10 N/m，回轉(zhuǎn)彈性系數(shù)kR=800 N/m，回轉(zhuǎn)阻尼系數(shù)cR=10 N/m，最大靜摩擦因子μ=0.08，時間步長=10-7s。此算例的目的是評估改進(jìn)模型處理固相行為的能力，實(shí)驗(yàn)中使用的玻璃珠密度為2 650 kg/m3。本實(shí)驗(yàn)中使用的粒徑為0.001 m，粒徑與實(shí)驗(yàn)所使用的玻璃微珠直徑保持一致，最終算例使用的粒子總數(shù)為4 108，算例是二維的實(shí)驗(yàn)1∶1建模，粒子數(shù)量由實(shí)驗(yàn)參數(shù)確定，這里不考慮粒子數(shù)量對模擬結(jié)果的影響。算例計(jì)算使用的時間步長均在符合Courant-Friedrichs-Lewy條件的前提下，考慮計(jì)算資源的限制得到。由于玻璃微珠不是流體，黏性系數(shù)無法直接得到，本文中所使用的數(shù)值0.08是參考摩擦系數(shù)得到的。

圖1 算例1計(jì)算域初始粒子排布

圖2示出算例1模擬結(jié)果，包括粒子配置，并與對應(yīng)的無量綱時間點(diǎn)的實(shí)驗(yàn)快照進(jìn)行了比較，無量綱時間T通過下式計(jì)算。

(18)

式中：t為真實(shí)時間；h0為固體滑坡在t=0時底端和頂端的距離。

由圖2可見，實(shí)驗(yàn)與數(shù)值模擬結(jié)果顯示出相似的特征，總體符合得較好。對于圖2a中存在的模擬結(jié)果與實(shí)驗(yàn)結(jié)果符合不是太好的問題，是由于實(shí)驗(yàn)開始時，擋板并不是瞬時移出，滑坡會受到擋板的干擾。對于圖2b中出現(xiàn)的模擬粒子積分面積小于實(shí)驗(yàn)顆粒面積的情況，是因?yàn)閷?shí)驗(yàn)是三維的，固體滑坡在拍攝方向上的剖視圖面積不均勻，而拍攝得到的面積必然會大于平均面積，也就顯得積分面積大于二維模擬結(jié)果的積分面積。顆粒材料的運(yùn)動也通過繪制滑坡前沿?zé)o量綱滑動的長度來量化和驗(yàn)證，如圖3所示，實(shí)驗(yàn)測得的滑坡前沿?zé)o量綱滑動長度與模擬得到的結(jié)果進(jìn)行對比，移動前沿的相對誤差在固體顆粒運(yùn)動穩(wěn)定后，最大不超過5.5%。

圖2 算例1模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比

圖3 算例1滑坡前沿?zé)o量綱滑動的長度對比

2.2 水下滑坡算例

算例2計(jì)算域初始粒子設(shè)置如圖4所示，這是一個水下滑坡的三維模型，在移除擋板后玻璃粉末開始運(yùn)動。算例2的幾何參數(shù)為：z=0.12 m，h1=0.4 m，h2=0.334 m，h3=0.32 m，α=35°。DEM計(jì)算中需要使用的參數(shù)為：kE=105N/m，cE=50 N/m，kN=104N/m，cN=10 N/m，kR=800 N/m，cR=10 N/m，μ=0.08，時間步長=10-5s。算例2基于文獻(xiàn)[20]中浸沒在水中的顆粒材料(玻璃粉末)沿35°光滑斜面滑動實(shí)驗(yàn)設(shè)計(jì)，算例2能評估模型在模擬大量離散固體流固耦合問題的能力。

圖4 算例2計(jì)算域初始粒子排布

將所建立的流固耦合模型和初始模型數(shù)值模擬結(jié)果與文獻(xiàn)[19]中的實(shí)驗(yàn)結(jié)果進(jìn)行比較，如圖5所示。實(shí)驗(yàn)結(jié)果用高速攝影機(jī)記錄，記錄時刻為0.02、0.17、0.32和0.47 s。在0.17 s時，流體表面的結(jié)構(gòu)與計(jì)算結(jié)果相似。隨著玻璃粉末組成的滑坡開始下降，流體表面出現(xiàn)凹陷，這在計(jì)算結(jié)果中也可清楚地觀察到。在0.32 s時，由于反彈，流體表面出現(xiàn)兩組波浪，這一現(xiàn)象在模擬計(jì)算結(jié)果中也得到了清晰的反映。在0.47 s時，剖面圖顯示此剖面上滑坡并不連續(xù)，但整個三維的滑坡是連續(xù)的，即在z方向的其他位置是連續(xù)的，這是由于MPS方法模擬的固體粒子尺寸與實(shí)際玻璃粉末尺寸相比較大，使得固體粒子數(shù)相對較少，隨著滑坡的進(jìn)行，部分位置就會出現(xiàn)這種情況，是合理的。由圖5可見，滑坡形狀在各時間點(diǎn)與實(shí)驗(yàn)結(jié)果對應(yīng)得較好，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，表明該模型具有分析流固耦合問題的能力。

根據(jù)滑坡移動前沿在斜坡上相對于初始前沿點(diǎn)的移動距離，定量分析了滑坡前緣的位置，結(jié)果如圖6所示，模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。在初始階段，由于擋板的移出對實(shí)驗(yàn)結(jié)果影響較大，距離變化較小，相對誤差較大，不具有參考價值。運(yùn)動穩(wěn)定后最大相對誤差僅為8.6%，證明了模型的準(zhǔn)確性。

3 結(jié)論

本文在基于MPS方法的核反應(yīng)堆關(guān)鍵熱工安全現(xiàn)象分析軟件平臺(PANDA)上開發(fā)了適用于分析流固耦合問題的PANDA-DEM模塊，并通過引入多相流的思路來處理相界面上的黏性來提高所使用模型的穩(wěn)定性和準(zhǔn)確性。對固體滑坡和水下滑坡算例進(jìn)行了數(shù)值模擬，并與實(shí)驗(yàn)結(jié)果進(jìn)行對比。通過定性對比固相的整體運(yùn)動和液相流場，定量比較滑坡前沿的移動距離，對模型進(jìn)行了驗(yàn)證，結(jié)果表明：1)本研究將離散單元法引入MPS方法框架內(nèi)構(gòu)建的DEM模型可用于模擬固體顆粒滑坡行為，滑坡前沿?zé)o量綱移動長度實(shí)驗(yàn)和模擬的相對誤差不超過6.3%；2)本研究中基于DEM模型和MPS方法構(gòu)建的流固耦合方法可用于模擬水下固體滑坡行為，滑坡前沿移動距離實(shí)驗(yàn)和模擬的相對誤差不超過8.6%。因此，本研究中建立的模型可對包含大量離散固體流固耦合問題進(jìn)行準(zhǔn)確分析，后續(xù)工作會將此模型擴(kuò)展到對核反應(yīng)堆堆芯熔化嚴(yán)重事故的分析模擬中。