基于MPS方法的離散固體與流體相互作用數值模擬研究

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

(1.西安交通大學 動力工程多相流國家重點實驗室，陜西 西安 710049；2.西安交通大學 核科學與技術學院，陜西 西安 710049)

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

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

在核反應堆中發生堆芯熔化嚴重事故時，堆芯熔化后產生的熔融物發生遷徙，未發生熔化的燃料微粒、基體碎片、包殼碎片也被熔融物夾帶著發生遷徙，這是一個包含大量離散固體的流固耦合問題。已有的嚴重事故程序對堆芯熔化問題進行研究時，大都采用集總參數的燭化模型，如MELCOR等，這種模型計算得到的結果較為粗糙，且不考慮未熔化的離散固體之間的相互作用。

本研究使用MPS方法對這種現象進行模擬，對熔融物夾帶的固體之間的相互作用以及液相與固相之間的作用進行分析。現有的處理方法難以進行準確的模擬，需要進行相關模型的開發。

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

1 數值模型

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

1.1 液相

使用MPS方法來對流體的運動進行模擬，控制方程包括連續性方程和動量守恒方程：

(1)

(2)

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

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

1.2 固相

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

DEM的主要思路是對每一個離散單元進行受力分析，通過牛頓第二定律得到其運動狀態。

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

(3)

(4)

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

(5)

(6)

(7)

(8)

式中，ΔtDEM為DEM計算的時間步長。

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

(9)

(10)

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

(11)

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

(12)

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

1.3 流固耦合

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

(13)

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

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

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

(14)

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

(15)

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

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

(16)

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

此時黏性項采用式(17)進行計算：

G(‖rj-ri‖,re)

(17)

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

一般情況下，MPS方法計算的時間步長要大于DEM計算，時間步長的耦合通過先對流體粒子推進一個MPS方法的時間步長，如果存在DEM計算，就進入DEM計算內循環，推進數個DEM時間步長，令其總的時間步長等于MPS方法計算時間步長來實現。

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

2 滑坡算例的模擬

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

2.1 固體滑坡算例

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

圖1 算例1計算域初始粒子排布

圖2示出算例1模擬結果，包括粒子配置，并與對應的無量綱時間點的實驗快照進行了比較，無量綱時間T通過下式計算。

(18)

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

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

圖2 算例1模擬結果與實驗結果的對比

圖3 算例1滑坡前沿無量綱滑動的長度對比

2.2 水下滑坡算例

算例2計算域初始粒子設置如圖4所示，這是一個水下滑坡的三維模型，在移除擋板后玻璃粉末開始運動。算例2的幾何參數為：z=0.12 m，h1=0.4 m，h2=0.334 m，h3=0.32 m，α=35°。DEM計算中需要使用的參數為：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基于文獻[20]中浸沒在水中的顆粒材料(玻璃粉末)沿35°光滑斜面滑動實驗設計，算例2能評估模型在模擬大量離散固體流固耦合問題的能力。

圖4 算例2計算域初始粒子排布

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

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

3 結論

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