嚴重事故熔池關(guān)鍵現(xiàn)象的MPS模型優(yōu)化

朱影子,熊進標

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

核反應堆發(fā)生失水事故后,堆內(nèi)冷卻劑大量喪失,導致燃料元件冷卻不足,發(fā)生堆芯熔化事故[1-3]。當堆芯熔化過程發(fā)展到一定程度,堆芯熔融物將落入壓力容器的下腔室,在瞬態(tài)熔池形成過程中,未熔化的堆芯材料,包括:鋯合金包殼,不銹鋼結(jié)構(gòu)件,銦科鎳定位格架等在燃料衰變熱的加熱下繼續(xù)熔化。熔化后的材料間存在密度差,則會發(fā)生密度分層現(xiàn)象,進而影響壓力容器下封頭熱流密度的分布,對其結(jié)構(gòu)完整性造成威脅。因此,有必要探索熔池內(nèi)的固體塊熔化分層的機理,為核反應堆安全分析提供相應的理論支持。

關(guān)于固體的熔化現(xiàn)象,已經(jīng)有很多的研究,其中Guo等[4]運用有限體積粒子法模擬了伍德合金的熔化,并與實驗作對比,模擬結(jié)果與實驗吻合較好。Li 等[5]對兩種密度不同、粘性不同的不相容液體的分層做了實驗和模擬,二者結(jié)果較為吻合。在熔化分層現(xiàn)象的模擬中,固液相界面發(fā)生較大的變化,因此選用Koshizuka[6]提出的基于拉格朗日描述的移動粒子半隱式法(MPS),此方法沒有網(wǎng)格約束,在計算域使用自由移動的粒子,便于處理大變形流動問題。此方法在核反應堆中[7]已應用于傳熱相變[8,9],液滴沖擊[10],熔融物遷移[11]和熔融物-混凝土相互作用[12]等問題的模擬。

熔化分層現(xiàn)象涉及熔化、分層等多個物理過程,在熔化過程中使用基于焓值相變模型來模擬,并且用黏性模型來計算固、液之間的黏性變化,其中改進黏性顯式算法為隱式算法來保證計算的高效性和穩(wěn)定性,運用PMS模型來實現(xiàn)固液耦合,并將熔化模擬結(jié)果和Guo[4]的實驗作對比,二者結(jié)果吻合較好。同時對兩種不同密度、不同黏性的流體做了分層模擬,模擬結(jié)果與Li[5]的實驗對比,結(jié)果顯示較好的一致性。最終將熔化現(xiàn)象和分層現(xiàn)象結(jié)合起來,模擬了伍德合金固體在密度較大的液體中熔化并分層的現(xiàn)象,并預測其熔化分層過程。

1 MPS方法

熔化過程中的控制方程包括連續(xù)性方程、動量方程和能量方程:

式中:ρ、ν、k——密度、運動黏度、導熱系數(shù)。

方程(2)中右側(cè)第三項為固體粒子和液體粒子之間的相互作用力,可以通過壓力差近似得到:

固體的控制方程為:

在粒子法中,控制方程的離散是通過粒子間相互作用模型實現(xiàn)的,如圖1所示。其中核函數(shù)是最基本的模型,如公式(7)所示,是用來計算粒子數(shù)密度的一個權(quán)重函數(shù),有很多類型,這里采用Koshizuka[6]提出的模型,如圖2。

式中:re——有效半徑,只有作用域內(nèi)的粒子之間有相互作用,此值一般取粒子直徑的2～3倍。

圖1 粒子有效半徑Fig.1 The cut-off radius

圖2 核函數(shù)Fig.2 Kernel function

梯度算子基于粒子間梯度權(quán)重平均由梯度模型離散,如公式(8)及圖3所示。

式中:N d——空間維數(shù);

n0——初始標準粒子數(shù)密度。

為了計算的穩(wěn)定性,壓力計算采用了Kondo[13]改進后的梯度模型,如公式(9)所示。

式中:β、γ——與時間步長、粒子大小等有關(guān)的修正系數(shù)。

β和γ應滿足公式(10)的邊界條件范圍

圖3 梯度模型Fig.3 The gradient model

拉普拉斯模型的物理示意如圖4所示。

圖4 拉普拉斯模型Fig.4 The Laplacian model

其中λ為修正計算結(jié)果,如公式(11)所示。

2 熔化模型

熔化是一個涉及相變的復雜物理過程,本文使用基于焓值的熔化模型來計算固-液相變過程,為了計算穩(wěn)定,使用黏性模型計算固、液相變之間的黏性,應用PMS模型實現(xiàn)固液之間的耦合。

2.1 基于焓的熔化相變模型

固體熔化過程中的相變是由公式(3)計算,溫度T和液相份額ε與焓值H的關(guān)系如公式(12)～公式(13)及圖5所示。

式中:C p,s和C p,l——熔化前后固相和液相的比熱。

圖5 基于焓值的相變模型Fig.5 The phase change model based on enthalpy

2.2 黏性模型

在實際熔化過程,固體的黏性為無窮大,遠遠大于液體,為了計算的穩(wěn)定性,引入Ramacciotti[14]的黏性模型。

式中:νi——運動黏性;

Hεcr——固體粒子變?yōu)橐后w粒子的臨界焓值,一般取εcr為0.375時的焓值;

A——流變 系 數(shù),根 據(jù)Guo[4]的 研究,A取值為0.054。

固體的黏性較大,要滿足擴散穩(wěn)定性,計算時間步長非常小,為了提高計算效率,將黏性項的顯式計算改為隱式計算:

2.3 流固耦合模型

熔化過程中涉及固體向液體的轉(zhuǎn)變,因此要處理流固耦合的問題,使用PMS模型[15]可以模擬此現(xiàn)象。其計算理論是:先將所有粒子視為液體粒子通過公式(1)～公式(3)計算其速度和位置,解能量方程由焓值判定識別固體粒子,計算固體塊的速度和重心位置,最后再對固體粒子的速度和位置做修正(見公式(17)～公式(22))。

其中N是粒子總數(shù)。

粒子的幾何重心為:

粒子的角速度:

每個固體粒子具有相同的密度,轉(zhuǎn)動慣量I為:

最后,修正固體粒子的速度和位置:

最終,計算的流程圖如圖6所示。

圖6 計算流程圖Fig.6 The workflow of the calculation

3 結(jié)果

3.1 熔化模擬

為了驗證基于焓值的熔化模型的準確性,模擬了伍德合金的熔化,并與Guo[4]的實驗作對比,將伍德合金置于恒定溫度的加熱板上,并監(jiān)測熔化過程中固體塊的高度和液膜擴展長度的變化。伍德合金的物理性質(zhì)如表1所示,計算中忽略固、液之間密度的變化。

表1 伍德合金物性Table 1 Properties of the Wood alloy

續(xù)表

初始幾何設(shè)置如圖7所示,開始所有粒子均為固體粒子,初始溫度均勻分布,為27℃,進行了底部加熱板的溫度分別為250℃,300℃和350℃的模擬。其他物理量設(shè)置均與實驗相同。熔化過程持續(xù)2 s,以導熱為主,此熔化過程持續(xù)時間較短,計算中忽略與空氣的對流換熱和輻射換熱。模擬結(jié)果如圖8和圖9所示,與實驗值相比固體塊的高度變化和液膜的擴展均偏低,除了計算本身的誤差以外,主要為計算中忽略輻射和對流換熱。

圖7 熔化模擬幾何Fig.7 Geometry of the Wood alloy

圖8 液膜擴展長度的變化Fig.8 Variation of the spreading width

圖9 固體高度的變化Fig.9 Variation of cube height

3.2 分層模擬

由于密度差產(chǎn)生分層是熔池形成行為中重要現(xiàn)象。Li[5]針對兩種互不相溶的液體做了分層實驗。實驗描述如圖10所示,初始時刻,中間有一塊隔板將兩種液體隔開,實驗開始后,隔板以0.2 m/s 的速度向上運動。鹽水的密度和黏性由其濃度決定,兩種液體界面處的密度和黏性由二者的算數(shù)平均值獲得。

圖10 分層模擬幾何Fig.10 Geometry for stratification

實驗分為兩組,A組使用黏性較大的硅油(5000 mm2/s)和10%(質(zhì)量分數(shù))的NaCl水溶液,B組使用黏性較小的硅油(200 mm2/s)和20%(質(zhì)量分數(shù))的NaCl水溶液。實驗中為了區(qū)分二者,在NaCl水溶液中添加了亞甲藍染色劑,具體參數(shù)如表2所示。

表2 兩種分層液體的參數(shù)Table 2 Properties of the fluids in stratification

圖11 分層實驗與模擬結(jié)果對比Fig.11 Comparison of the simulation and experiment

模擬結(jié)果如圖11(a)顯示了黏性比較大的硅油和鹽水的分層實驗與模擬的對比結(jié)果,圖中密度較小的硅油上浮,密度較大的鹽水下沉,并最終產(chǎn)生穩(wěn)定分層的結(jié)果。在前5 s,液體流動較快,且在容器底部出現(xiàn)一層硅油層,之后液體流動較為緩慢,分層形狀基本不變。這是因為隔板初始抽出之后,兩種液體界面有較大的壓力差驅(qū)動產(chǎn)生較大的動能,此能量隨著時間而衰減。圖11(b)顯示了黏性比較小的兩種液體的分層結(jié)果,與A組相比,B組實驗中硅油的黏性較小,鹽水的密度較大,當隔板抽出后,兩種液體的界面處產(chǎn)生瞬時震蕩現(xiàn)象,模擬中也將這種趨勢預測出。

3.3 熔化分層模擬

將熔化現(xiàn)象和分層現(xiàn)象結(jié)合起來,模擬伍德合金固體在密度較大的液體(稱為背景液體)中熔化并分層的現(xiàn)象。其中伍德合金固體密度為8 528 kg/m3,并且固體熔化前后密度保持不變,背景液體密度為9 528 kg/m3。初始幾何如圖12所示,設(shè)計固體在背景液體中的位置處于其所受浮力和其重力剛好抵消,固體和液體的相對速度較小,計算中簡化固體在液體中的熔化現(xiàn)象,忽略固、液之間的對流換熱,以導熱為主。假設(shè)伍德合金熔化后的液體和背景液體有相同的比熱、導熱系數(shù)、粘性等物理參數(shù)。初始條件:背景液體溫度為350℃,固體溫度為27℃;邊界條件:幾何左右兩側(cè)為絕熱邊界,底部為定溫350℃。伍德合金固體熔化后的液體和背景液體的界面處,密度為二者的平均值。

圖12 熔化分層模擬初始幾何Fig.12 Geometry for melting and stratification

伍德合金在背景溶液中熔化分層結(jié)果的預測如圖13所示,固體在63.876 s完全熔化,并穩(wěn)定分層。在0～60 s中,使用不同相態(tài)區(qū)分粒子,觀察固體形狀的變化,淺藍色為固體粒子,深藍色為液體粒子,固體被加熱,并逐漸熔化。65 s時,固體完全熔化并穩(wěn)定分層,使用不同物質(zhì)區(qū)分粒子,其中上層紅色為密度較小的伍德合金熔化液體,下層藍色為密度較大的背景液體。最終液體的壓力分布符合靜壓分布規(guī)律,如圖14所示。圖15為最終的溫度分布,固體熔化后的液體處溫度最低,以熔化后的液體為中心,周圍溫度逐漸升高,底部邊界維持350℃的恒溫。模擬結(jié)果符合熔化分層的物理規(guī)律。

圖13 熔化分層模擬結(jié)果Fig.13 Process of melting and stratification

圖14 在65 s時的壓力分布Fig.14 Pressure distribution at 65 s

圖15 在65 s時的溫度分布Fig.15 Temperature distribution at 65 s

4 結(jié)論

為了探索核反應堆嚴重事故后,在瞬態(tài)熔池形成過程中,燃料元件等在熔池中的熔化及分層現(xiàn)象。將原始的移動粒子半隱式法(MPS)進行改進,首先將壓力源項改為混合源項提高壓力計算的穩(wěn)定性,并將黏性顯式計算改為隱式算法,提高計算的效率。然后加入能量方程并對熔化過程進行建模,結(jié)合基于焓值的熔化/凝固相變模型、黏性模型、流固耦合模型等模擬了三維伍德合金的熔化現(xiàn)象,模擬結(jié)果與實驗結(jié)果較為吻合。通過模擬互不相溶的硅油和鹽水驗證算法對分層現(xiàn)象的有效性,并與實驗結(jié)果進行對比,顯示較好的吻合度。最終將熔化和分層現(xiàn)象結(jié)合,模擬了伍德合金在背景液體中的熔化分層現(xiàn)象,并預測其過程,通過對結(jié)果的分析,表明該模擬符合熔化分層的物理規(guī)律。