朱光昱， 靖劍平， 石興偉， 左嘉旭， 劉宇生，4， 溫爽*

(1.生態(tài)環(huán)境部核與輻射安全中心， 北京 100082； 2.國(guó)家環(huán)境保護(hù)核與輻射安全審評(píng)模擬分析與驗(yàn)證重點(diǎn)實(shí)驗(yàn)室， 北京 102488; 3.中國(guó)核電工程有限公司， 北京 100084; 4.哈爾濱工程大學(xué)核科學(xué)與技術(shù)學(xué)院， 哈爾濱 150001)

核電廠嚴(yán)重事故中，堆芯構(gòu)件由于失去冷卻而熔化形成高溫熔融池[1]。在中國(guó)的第三代先進(jìn)壓水堆發(fā)展過(guò)程中，傾向于采用堆內(nèi)熔融物滯留技術(shù)(in-vessel retention, IVR)作為核電廠嚴(yán)重事故的緩解措施[2]，該策略通過(guò)對(duì)反應(yīng)堆壓力容器(reactor pressure vessel, RPV)外壁面進(jìn)行冷卻帶走堆內(nèi)熔融池產(chǎn)生的衰變熱，從而保持 RPV 的結(jié)構(gòu)完整性。在高溫條件下，RPV下封頭會(huì)發(fā)生明顯的熱膨脹，從而導(dǎo)致其外部冷卻流道結(jié)構(gòu)與初始設(shè)計(jì)存在偏差[3]。如果上述變化導(dǎo)致局部傳熱惡化，熔融池導(dǎo)出熱量大于局部CHF值，則可能發(fā)生RPV熔穿進(jìn)而導(dǎo)致IVR策略失效，因此在IVR策略設(shè)計(jì)過(guò)程中，需要對(duì)嚴(yán)重事故后RPV的形變進(jìn)行研究。

在高溫熔融池的作用下，RPV壁面會(huì)發(fā)生熔化進(jìn)而降低RPV的壁面厚度，最終加劇了RPV 壁面的變形。因此在對(duì)RPV下封頭的熱膨脹變形影響進(jìn)行評(píng)估時(shí)，也必須考慮鋼結(jié)構(gòu)熔化帶來(lái)的影響。在熔融池的演化過(guò)程中，由于其中的輕金屬與氧化物不相溶，最終會(huì)形成雙層熔融池結(jié)構(gòu)[4]。頂部較輕的金屬層導(dǎo)熱性能良好，產(chǎn)生的熱聚集效應(yīng)會(huì)使與金屬層接觸的RPV下封頭內(nèi)壁面的熔化可能更加明顯，在研究過(guò)程中同時(shí)需要考慮熔融池分層帶來(lái)的影響。

在工程設(shè)計(jì)上，一般通過(guò)數(shù)值模擬的方式對(duì)嚴(yán)重事故后RPV下封頭形變進(jìn)行分析。溫爽等[5]和田欣鷺等[6]基于ANSYS研究了熱膨脹和高溫蠕變對(duì)下封頭失效的影響。羅娟等[7]通過(guò)彈塑性模型分析了下封頭自重、熔融池內(nèi)壓和外部冷卻水壓力對(duì)下封頭變形的影響。劉兆陽(yáng)等[8]采用ABAQUS軟件對(duì)示范快堆熔融物收集裝置進(jìn)行了安全分析。高永建等[9]在完成溫度場(chǎng)分析的基礎(chǔ)上，通過(guò)Larson-Miller參數(shù)積累損傷理論對(duì)RPV蠕變進(jìn)行了研究。上述工作中均采用了將熔融池內(nèi)部與下封頭解耦的方式完成力學(xué)分析，可能引入一定的不確定性。對(duì)此，現(xiàn)采用多物理場(chǎng)耦合仿真平臺(tái)COMSOL，建立一個(gè)熱工水力和固體力學(xué)耦合計(jì)算模型，對(duì)單層和雙層熔融池結(jié)構(gòu)下RPV下封頭的形變進(jìn)行研究，為嚴(yán)重事故下RPV外壁面的換熱情況研究提供必要的輸入?yún)?shù)。

1 計(jì)算模型

1.1 計(jì)算域設(shè)置

圖1所示為參考HPR1000堆型設(shè)計(jì)參數(shù)建立的二維軸對(duì)稱雙層熔融池計(jì)算模型，熔融物區(qū)域設(shè)置為流體域，其中氧化物層高度為1.58 m，金屬層高度為0.56 m。下封頭區(qū)域設(shè)置為固體域，下封頭內(nèi)徑為4.4 m，下封頭壁面厚度為0.15 m。圖1中r表示軸對(duì)稱模型的半徑，g為重力加速度，模型中取9.8 m/s2。單層熔融池計(jì)算模型與圖1類似，但不存在分層結(jié)構(gòu)，總高度為1.8 m。為了減少邊緣效應(yīng)對(duì)力學(xué)計(jì)算的影響，下封頭頂部區(qū)域延伸了長(zhǎng)度為3 m的不銹鋼壁面，圖1中并未完全示出。

當(dāng)前COMSOL計(jì)算模型中，熔融池和下封頭的溫度場(chǎng)采用固體和流體傳熱模塊計(jì)算，熔融池湍流場(chǎng)采用低雷諾數(shù)k-ε模型計(jì)算，下封頭的熱膨脹過(guò)程采用固體力學(xué)模塊計(jì)算。熔融物和下封頭的相變過(guò)程采用相變模型計(jì)算。計(jì)算域中，金屬層頂部和下封頭延伸部分的內(nèi)壁面均設(shè)置為表面對(duì)環(huán)境輻射條件，發(fā)射率分別設(shè)置為0.45和0.8[4]。下封頭外側(cè)設(shè)置為400 K的恒溫壁面[5]。對(duì)于雙層熔融池結(jié)構(gòu)，可以認(rèn)為衰變熱僅存在于氧化物層中，計(jì)算公式[4]為

Q=0.125fdP0t-0.28

(1)

式(1)中：Q為堆芯衰變熱功率，MW；fd為衰變熱修正系數(shù)，設(shè)置為0.75[12]；P0可認(rèn)為是反應(yīng)堆正常運(yùn)行時(shí)的熱功率，本文取3 050 MW；t為停堆時(shí)間。

圖1 計(jì)算域和網(wǎng)格劃分示意Fig.1 Calculation domain and mesh generation

1.2 熔融池物性設(shè)置

熔融池內(nèi)通過(guò)引入相變模型模擬凝固過(guò)程，材料密度采用Boussinesq近似處理。參考VULCANO熔融物鋪展實(shí)驗(yàn)相關(guān)數(shù)值模擬經(jīng)驗(yàn)[10]，熔融池中糊狀區(qū)的黏度為

μ=μle2.5Cθs

(2)

式(2)中：μl為液相黏度；系數(shù)C一般取3.5～8[11]；θs為固相體積分?jǐn)?shù)，可直接在計(jì)算過(guò)程中實(shí)時(shí)調(diào)用參數(shù)。

單層和雙層熔融池內(nèi)其他計(jì)算所需物性參照文獻(xiàn)[12-14]設(shè)置，膨脹系數(shù)參考液態(tài)金屬和陶瓷氧化物設(shè)置，具體如表1和表2所示。

表1 單層熔融池物性參數(shù)表

表2 雙層熔融池物性參數(shù)表

1.3 壓力容器物性設(shè)置

計(jì)算過(guò)程中依然采用相變材料模型計(jì)算下封頭熔化過(guò)程，并為此設(shè)置了10 K的相變區(qū)間，RPV的具體物性設(shè)置如表3所示。

為了引入固體力學(xué)模型，下封頭需要設(shè)置為固體計(jì)算域，這導(dǎo)致無(wú)法直接計(jì)算RPV壁面熔化后混入熔融物的過(guò)程，因此對(duì)下封頭物性進(jìn)行了如下設(shè)置。通過(guò)預(yù)計(jì)算判斷，當(dāng)前熔融池參數(shù)設(shè)置下，僅下封頭與金屬層接觸區(qū)域會(huì)發(fā)生熔化，而熔化后的鋼液相對(duì)于金屬層總質(zhì)量而言很小，因此可以將下封頭熔化后區(qū)域的材料熱工物性參數(shù)設(shè)置為與金屬層一致。

此外，對(duì)于下封頭熔化區(qū)域，在材料本身的導(dǎo)熱系數(shù)之外將湍流導(dǎo)熱添加至熔化后區(qū)域的導(dǎo)熱系數(shù)當(dāng)中，湍流導(dǎo)熱的計(jì)算公式為

(3)

式(3)中：kT為湍流流動(dòng)產(chǎn)生的導(dǎo)熱系數(shù)；PrT為湍流普朗特?cái)?shù)，取0.9；μ為湍流動(dòng)力黏度；Cp為比熱容，兩項(xiàng)均參照熔化后的金屬層參數(shù)設(shè)置。

PRV的楊氏模量通過(guò)插值函數(shù)設(shè)置[7]，具體如表4所示。對(duì)于固體力學(xué)計(jì)算，將熔化區(qū)域的力學(xué)參數(shù)直接設(shè)置為0會(huì)導(dǎo)致計(jì)算發(fā)散。對(duì)此，對(duì)表3中熔化后的泊松比和表4中液相線以上的楊氏模量設(shè)置了極小的數(shù)值，以達(dá)到忽略熔化區(qū)域力學(xué)影響的目的。

表3 壓力容器物性參數(shù)表

表4 不同溫度下的楊氏模量

1.4 網(wǎng)格敏感性分析

采用瞬態(tài)方式求解，通過(guò)設(shè)置較長(zhǎng)的計(jì)算時(shí)長(zhǎng)獲得各物理場(chǎng)達(dá)到穩(wěn)定狀態(tài)時(shí)的參數(shù)?？紤]停堆3 h后已經(jīng)形成雙層熔融池結(jié)構(gòu)，此時(shí)氧化物層的衰變熱為21.23 MW。圖2所示為不同網(wǎng)格數(shù)計(jì)算得到的RPV下封頭外壁面應(yīng)力分布。當(dāng)網(wǎng)格加密至30 000以上后，網(wǎng)格變化對(duì)計(jì)算得到的應(yīng)力分布變化的影響降低至 1 %以內(nèi)，雙層熔融池后續(xù)計(jì)算采用的網(wǎng)格數(shù)為32 298。對(duì)于單層熔融池計(jì)算模型，采用相同的方法進(jìn)行網(wǎng)格敏感性分析，最終采用的網(wǎng)格數(shù)為21 279。

圖2 不同網(wǎng)格數(shù)下封頭外壁面應(yīng)力Fig.2 Stress distribution on the outer surface of lower head under under different mesh number

2 數(shù)值模擬結(jié)果分析

2.1 熔融池溫度場(chǎng)分析

圖3 單層熔融池和雙層熔融池的溫度分布 Fig.3 Temperature distribution of ingle-layer corium pool and two-layer corium pool

計(jì)算過(guò)程中，熔融池以參考溫度作為初始化溫度，下封頭的初始溫度設(shè)置為450 K。圖3為衰變熱為23.78 MW時(shí)的單層熔融池和衰變熱為21.23 MW時(shí)的雙層熔融池溫度分布情況。本文中，相位角指的是與計(jì)算域中對(duì)稱軸的夾角。對(duì)于單層熔融池結(jié)構(gòu)，下封頭中相位角大于25°的區(qū)域均不同程度的超過(guò)了不銹鋼的熔點(diǎn)。雙層熔融池中，與氧化物層接觸的下封頭區(qū)域溫度較低不會(huì)超過(guò)熔點(diǎn)，對(duì)于與金屬層區(qū)域接觸的區(qū)域，在熱聚集效應(yīng)作用下，這部分下封頭內(nèi)壁溫度會(huì)顯著升高并超過(guò)不銹鋼的熔點(diǎn)。

圖4 單層熔融池和雙層熔融池的下封頭von Mises應(yīng)力分布 Fig.4 Distribution of von Mises stress in lower head of single-layer corium pool and two-layer corium pool

2.2 RPV熱膨脹形變分析

圖4為單層熔融池和雙層熔融池狀態(tài)下封頭的von Mises應(yīng)力分布情況，高溫作用下不銹鋼壁面的熔化情況也可在圖中看出?？梢娤路忸^熱膨脹產(chǎn)生的應(yīng)力主要集中在壓力容器外壁面處。對(duì)于發(fā)生熔化的區(qū)域，在較高的溫度和壁面厚度由于熔化削弱的雙重作用下，下封頭外側(cè)von Mises應(yīng)力可達(dá)其他區(qū)域的兩倍以上。

圖5為單層熔融池和雙層熔融池狀態(tài)下RPV的熱膨脹形變量分布，為了便于觀察，圖5所示的膨脹量顯示為實(shí)際值的5倍。計(jì)算域中，RPV延伸區(qū)域頂部3 m處的邊界條件為縱向固定約束，導(dǎo)致計(jì)算域內(nèi)所有的縱向熱膨脹形變會(huì)逐漸向下累積。同時(shí)，由于對(duì)稱軸相當(dāng)于橫向位移約束，所以橫向形變的則會(huì)逐漸向下封頭的高相位角積累。在各處的局部熱膨脹形變和上述兩種積累效應(yīng)的共同作用下，單層熔融池的最大形變出現(xiàn)在相位角40°左右的區(qū)域，最大形變量達(dá)到19.7 mm。雙層熔融池的最大形變量出現(xiàn)在與金屬層接觸的區(qū)域，其中最大的膨脹位移達(dá)到了16.9 mm。對(duì)于下封頭的底部區(qū)域，雖然兩種熔融池狀態(tài)下的局部的熱膨脹形變均很小，但在縱向的積累效應(yīng)作用下依然存在約15 mm的膨脹位移。

2.3 RPV內(nèi)壓影響分析

除衰變熱和下封頭壁面熔化外，RPV內(nèi)壓也是其形變的重要影響因素[5-7]，在泄壓系統(tǒng)的作用下，可以認(rèn)為熔融池形成時(shí)，RPV的內(nèi)壓一般不會(huì)超過(guò)1 MPa。本文選取雙層熔融池結(jié)構(gòu)對(duì)內(nèi)壓影響進(jìn)行分析，通過(guò)施加邊界載荷改變內(nèi)壓。圖6所示為不同衰變熱和內(nèi)壓條件下的下封頭外壁面膨脹位移隨相位角的變化。圖6中衰變熱21.23 MW和14.39 MW分別對(duì)應(yīng)停堆后約3 h和12 h的時(shí)刻。衰變熱為21.23 MW時(shí)，1 MPa和2 MPa內(nèi)壓分別使得最大熱膨脹位移增加了約0.9 mm和1.4 mm，這將導(dǎo)致RPV對(duì)冷卻流道的擠壓更加明顯，使外流道傳熱性能發(fā)生變化，有可能造成冷卻惡化和IVR策略的失效。隨著衰變熱降低，下封頭各相位角的徑向膨脹位移呈現(xiàn)降低的趨勢(shì)，然而當(dāng)內(nèi)壓較高時(shí)這種收縮較為緩慢，從而導(dǎo)致RPV冷卻流道變形時(shí)間增長(zhǎng)。

圖5 單層熔融池和雙層熔融池的RPV的熱膨脹形變量 Fig.5 RPV thermal expansion deformation of single-layer corium pool and two-layer corium pool

圖6 下封頭外側(cè)徑向熱膨脹位移Fig.6 Thermal expansion deformation displacement of lower head outboard

3 結(jié)論

根據(jù)HPR1000堆型設(shè)計(jì)參數(shù)，采用COMSOL軟件搭建了一個(gè)耦合計(jì)算模型，對(duì)嚴(yán)重事故后單層和雙層熔融池結(jié)構(gòu)下的RPV熱膨脹形變進(jìn)行了分析。得出以下結(jié)論。

(1)嚴(yán)重事故后RPV內(nèi)的熔融池將導(dǎo)致IVR下封頭內(nèi)壁面局部區(qū)域發(fā)生熔化，受局部高溫和不銹鋼壁面變薄的共同影響，RPV外冷卻流道受到明顯擠壓。對(duì)于單層熔融池結(jié)構(gòu)，最大形變出現(xiàn)在相位角40°左右的區(qū)域，最大形變量達(dá)到19.7 mm。對(duì)于雙層熔融池結(jié)構(gòu)，最大形變出現(xiàn)在相位角80°左右的區(qū)域，最大形變量達(dá)到16.9 mm。

(2)RPV膨脹形變量隨著內(nèi)壓增高而明顯增加，由此可能導(dǎo)致外流道傳熱性能發(fā)生變化，進(jìn)而造成冷卻惡化和IVR策略的失效，在工程設(shè)計(jì)過(guò)程中應(yīng)予以關(guān)注。