林 萍,楊森皓,陳旭鵬,銀建中,韓志遠,謝國山

(1.大連理工大學(xué) 化工學(xué)院,遼寧大連 116024;2.中國特種設(shè)備檢測研究院,北京 100029)

0 引言

安全性是核電研究的重中之重,反應(yīng)堆壓力容器(RPV)作為核電站中包容反應(yīng)堆活性區(qū)的重要設(shè)備和安全屏障,長期在高溫高壓高輻照的惡劣工況下服役。當(dāng)失水事故(LOCA)發(fā)生時,反應(yīng)堆快速停堆,主冷卻劑流量降低,應(yīng)急堆芯冷卻系統(tǒng)(ECCS)啟動,低溫安注水大量注入環(huán)腔,產(chǎn)生較大溫度梯度,進而產(chǎn)生熱應(yīng)力,與壓力的共同作用下產(chǎn)生承壓熱沖擊(PTS)現(xiàn)象[1],使RPV內(nèi)壁產(chǎn)生裂紋缺陷并持續(xù)擴展或使原有缺陷出現(xiàn)擴展[2],甚至有泄漏風(fēng)險。因而如何更精確地評估PTS現(xiàn)象下RPV的完整性,是目前亟待解決的問題。

PTS工況下,反應(yīng)堆停堆后由于堆芯內(nèi)剩余裂變和中子俘獲產(chǎn)物衰變產(chǎn)生功率,其功率是一個逐漸下降的過程,開始時下降速度很快,達到一定數(shù)值后,下降速度變得緩慢。雖然此值只有穩(wěn)態(tài)功率的百分之幾,但仍然對壁面溫度的變化有阻礙作用[3-4],直接影響PTS分析結(jié)果。馬梓淇等[5]在對AP1000型號的RPV在PTS事故下的熱應(yīng)力數(shù)值模擬時,考慮了1 MW/m3堆芯熱功率的影響。楊森皓等[6]基于三回路RPV模型,進行了PTS瞬態(tài)模擬與斷裂分析,也設(shè)置了1 MW/m3的堆芯衰變熱,并對比了傳統(tǒng)無堆芯衰變熱的計算結(jié)果,得到兩種PTS事故工況下的衰變熱的影響情況。文獻[5-6]中均考慮了衰變熱的因素,但是所用的衰變熱絕對值偏小,且未考慮衰變熱是動態(tài)變化的過程。

在運行過程中,堆芯區(qū)域由于快中子長期輻照導(dǎo)致RPV材料脆化,韌脆轉(zhuǎn)變溫度升高[7]。為防止母材的腐蝕和脆化,在堆芯區(qū)域RPV內(nèi)壁有一層厚度3～10 mm的奧氏體不銹鋼材料作為堆焊層[8]。堆焊層的相關(guān)材料參數(shù)與母材材料參數(shù)不同,對分析結(jié)果有一定影響。

本文在前人研究基礎(chǔ)上,基于ACP1000三回路壓水堆,對25 cm2小破口失水事故(SBLOCA)工況進行計算;在熱工水力計算中,考慮停堆后變化的堆芯衰變熱的作用,分析事故下冷卻劑流動軌跡、RPV內(nèi)壁冷卻路徑,在冷點處假設(shè)裂紋進行斷裂力學(xué)分析;再通過與1 MW/m3堆芯衰變熱功率計算結(jié)果比較,得到變化堆芯衰變熱對同一工況的影響結(jié)果。

1 數(shù)值計算模型

1.1 物理模型

基于ACP1000三回路壓水堆,應(yīng)用三維繪圖軟件進行幾何建模。外部殼體設(shè)計高度為8 000 mm,外徑3 200 mm,壁厚200 mm,堆焊層厚度4 mm,殼體內(nèi)包裹堆芯區(qū)域及其他附件,由于其實際結(jié)構(gòu)較復(fù)雜,建模過程中進行了相應(yīng)的簡化。如圖1所示,其上方交錯設(shè)置3個進水管和3個出水管,其下方設(shè)置2個安注接管(DVI),DVI采用文丘里管形式。低溫安注水經(jīng)由DVI注入RPV環(huán)腔內(nèi),與腔內(nèi)原有流體混合,流經(jīng)下腔室,通過流量分配裙分配流量后,經(jīng)由堆芯支撐鍵流入堆芯、帶走熱量,再從熱段出水口流出。將幾何模型導(dǎo)入Design Modeler中,進行Slice劃分成多個比較規(guī)則的形狀,有利于后續(xù)劃分高質(zhì)量網(wǎng)格,并進行流體填充。進行網(wǎng)格劃分時,綜合考慮計算精度與計算要求,對DVI等處進行了網(wǎng)格加密處理,最終劃分網(wǎng)格數(shù)約300萬。

圖1 RPV物理模型Fig.1 Physical model of RPV

1.2 堆芯衰變熱計算

反應(yīng)堆停堆后功率包括剩余裂變產(chǎn)生的功率和中子俘獲產(chǎn)物產(chǎn)生的功率。停堆后的功率N(t)與停堆前功率N(0)相對變化可表示為式(1),式中右側(cè)第一項是剩余裂變功率,第二項是衰變功率[3]。

(1)

式中,t為停堆時間,s;n(0),n(t)為停堆前和停堆后t時刻的中子密度,中子/cm3;NS(t)為停堆后t時刻的衰變功率,W。

以低濃縮235U作為燃料的輕水堆相對中子密度隨時間的變化情況為:

(2)

本文選用的RPV正常運行功率3 180 MW,計算時假設(shè)停堆前反應(yīng)堆已經(jīng)運行了無限長的時間,由于在計算中忽略了其他俘獲產(chǎn)物的衰變功率,故將最終計算結(jié)果乘以1.1的安全系數(shù)。本文堆芯體積為8.375 64 m3,最終得到2 000 s內(nèi)單位體積的堆芯衰變熱功率如圖2所示。

圖2 2 000 s內(nèi)單位體積堆芯衰變熱功率變化情況Fig.2 Core decay heat power variation per unit volume in 2 000 s

1.3 控制方程

環(huán)腔內(nèi)流體假設(shè)為粘性不可壓縮湍流模型,其基本方程為Reynolds時均方程,選擇Realizablek-ε兩方程模型[9],此模型Cμ是與ε有關(guān)的公式。耗散率傳輸方程相較標(biāo)準(zhǔn)模型在生成項處不包含Pk,減少項中不具有奇異性。Realizablek-ε湍流模型控制方程[10]如下:

(3)

式中,μt為流體粘性系數(shù),kg/(m·s),μt=ρCμk2/ε,其中Cμ=(A0+ASUk/ε)-1;Pk,Gb為由平均速度梯度、浮力引起的湍流動能,J;σk,σε為湍流動能k、湍流耗散率ε的普朗特數(shù),σk=1.0,σε=1.2;C1,C2,C1ε為常數(shù),C1=max[0.43,η/(η+5)],其中η=Sk/ε。

流體域計算的溫度、壓力等參數(shù)與壁面的耦合邊界上,要滿足溫度和熱流密度的連續(xù)性條件:

(4)

式中,Tfp,Tsp為交界面兩側(cè)流體、固體的網(wǎng)格節(jié)點溫度;qf,qs為流體、固體兩側(cè)的熱流密度。

冷卻劑與RPV接觸,熱量傳遞主要形式為對流換熱及熱傳導(dǎo),對流換熱邊界條件為:

(5)

式中,h為流體和固體間的表面對流換熱系數(shù),W/(m2·℃);TS,Tf為RPV壁面溫度、流體溫度,℃。

固體傳熱方程為:

(6)

式中,ρ為材料密度,kg/m3;Cp為材料的定壓比熱容,J/(kg·℃);λ為材料的導(dǎo)熱系數(shù),W/(m·℃);qv為堆芯單位體積釋熱率,與堆芯熱功率有關(guān),W/m3。

1.4 求解設(shè)置

1.4.1 事故工況

本文計算的PTS工況是25 cm2小破口失水事故工況,該工況下的輸入溫度、壓力條件見圖3,在冷段進水口設(shè)置流量、溫度入口,熱段作為壓力出口。DVI接口安注水溫度為10 ℃,小破口失水事故中ECCS以壓力作為觸發(fā)條件,t=6.5 s觸發(fā)堆芯補水箱(CMT),在t=106 s時內(nèi)壓下降到蓄壓安注箱(ACC)的觸發(fā)條件,具體參數(shù)見表1。

圖3 25 cm2小破口失水事故輸入條件Fig.3 Input conditions for 25 cm2 small break loss of coolant accident

表1 冷卻劑注入各項參數(shù)Tab.1 Coolant injection parameters

1.4.2 條件設(shè)置

基于CFX平臺進行流體仿真計算,堆芯簡化為多孔模型[11-12],將冷卻劑、堆芯、RPV外殼、堆焊層劃分為4個區(qū)域,通過設(shè)置Fluid-Porous型、Fluid-Solid型和Solid-Solid型的交界面,實現(xiàn)不同區(qū)域之間的耦合。

RPV外殼和堆焊層兩者均設(shè)置為固體域,RPV殼體、接管和主要結(jié)構(gòu)材料為16MnNiMo(16MND5),堆焊層材料為奧氏體不銹鋼[13],材料物性參數(shù)分別如表2、表3所示。RPV外壁絕熱,材料計算模型應(yīng)用雙線性等向強化的彈塑性模型。冷卻劑為流體域,介質(zhì)為水,水的物性采用IAPWS-IF97標(biāo)準(zhǔn)計算,初始采用溫度300 ℃、壓力15.267 MPa時的值,在計算中考慮浮力的影響,水的物性隨溫度、壓力變化。堆芯材料為二氧化鈾和鋯錫合金陶瓷芯塊,堆芯為多孔域,流動損失系數(shù)0.39 m-1,橫向損失系數(shù)10 m-1,設(shè)置的熱源項為隨時間變化的函數(shù)(見圖2)。

表2 16MND5材料參數(shù)Tab.2 Material parameters of 16MND5

表3 堆焊層材料參數(shù)Tab.3 Material parameters of surfacing layer

2 結(jié)果分析

2.1 流動混合分析

由于安注水在環(huán)腔內(nèi)與原有冷卻劑混合的同時伴隨著熱量的交換,因此,通過提取環(huán)腔內(nèi)冷卻劑溫度變化情況來分析安注水注入環(huán)腔后的流動軌跡。安注水注入環(huán)腔后的溫度變化情況如圖4所示。

圖4 安注水在環(huán)腔內(nèi)流動混合情況Fig.4 The mixing condition of water flow in the annular cavity

當(dāng)只有CMT啟動時,安注水流速較小,安注外輪廓軌跡近似為一條直線;隨著安注時間變長,安注流體與環(huán)腔內(nèi)冷卻劑在環(huán)腔下部混合,且逐漸向整個環(huán)腔的周向區(qū)域擴散。t>106 s時,ACC啟動,安注流速增大,安注水流動軌跡近似為等腰三角形;隨著安注時間的增加,等腰三角形的頂角更大,環(huán)腔內(nèi)安注水大范圍擴散,環(huán)腔內(nèi)混合水溫開始逐漸降低,二者混合逐漸均勻。

安注水流經(jīng)環(huán)腔與冷卻劑混合后,經(jīng)由下腔室流入堆芯。從堆芯中部橫截面位置的溫度分布圖(見圖5(a))可以看出,堆芯內(nèi)水溫最高的區(qū)域集中在中心軸線附近,所以在堆芯中心線處沿A到B(從下到上,見圖5(b))提取軸向上的溫度。

圖5 提取堆芯溫度的路徑確定Fig.5 Determination of core temperature extraction path

圖6示出AB路徑上堆芯多孔域溫度的變化情況,對路徑長度進行了歸一化無量綱處理,即A點為0,B點為1。

圖6 冷卻劑沿堆芯軸線溫度變化Fig.6 Coolant temperature changes along the core axis

由于前50 s堆芯衰變熱功率較大,堆芯溫度不能完全排除,堆芯的加熱有一定的累加、滯后性。10 s時,低溫安注水僅有少量與原有冷卻劑混合后到達堆芯,堆芯主要冷卻劑還是溫度比較高的一回路冷卻水,此時堆芯熱功率大,水溫升高達到該壓力下的飽和溫度后,換熱能力下降,只有A點受冷卻較明顯。30 s時,伴隨著前期高熱功率的累積,堆芯溫度仍然在升高,但堆芯周圍冷卻劑逐漸循環(huán)帶走下方部分熱量,A,B點溫差最大為190.77 K左右,與理論流體吸收熱量計算公式Q=cmΔT所計算出的溫差相比誤差為6.41%,誤差在允許范圍內(nèi)[14]。時間遞進到110 s,ACC已經(jīng)啟動,大量低溫安注水注入RPV,堆芯冷卻速率加快,堆芯溫度較之前時刻降低,安注水流量增加,冷卻劑溫度降低,到達B點可換熱的冷卻液多,B點溫度下降。同條件下,計算1 MW/m3堆芯衰變熱作用下,B點相較A點溫度升高了30 K左右。

冷卻劑在流動過程中受壓力、浮力以及堆芯多孔域的影響,經(jīng)過堆芯的水流動軌跡存在一定的偏移現(xiàn)象,各處流體被流經(jīng)堆芯帶走熱量的路徑不可確定,所以不是越接近B點堆芯溫度越高。

安注冷卻劑在RPV環(huán)腔內(nèi)的流動情況會影響RPV內(nèi)壁的溫度分布,進而在內(nèi)壁產(chǎn)生較大熱應(yīng)力,所以在進行反應(yīng)堆壓力容器的完整性評估時,要根據(jù)安注冷卻劑流動軌跡,分析RPV內(nèi)壁上溫度的變化。

2.2 傳熱分析

ECCS啟動后,RPV內(nèi)壁面溫度場分布如圖7所示。低速安注初期,DVI下方有一條豎直的羽流區(qū)域[15-16],安注時間增大,羽流區(qū)域發(fā)生偏移,偏向DVI兩側(cè)。安注流量增大后,羽流區(qū)域更偏向于兩側(cè),此現(xiàn)象導(dǎo)致RPV環(huán)腔周向產(chǎn)生較大的溫度梯度,對環(huán)焊縫、缺陷強度和完整性評估有重要影響。

圖7 RPV內(nèi)壁溫度變化Fig.7 Inner wall temperature variation of RPV

進行完整性評估時,確定安注冷卻過程冷點的位置至關(guān)重要,冷點處溫度梯度大,產(chǎn)生的熱應(yīng)力也大,若此處存在缺陷,則最容易發(fā)生裂紋擴展和失效現(xiàn)象。為確定冷點位置,在RPV內(nèi)壁、堆焊層與母材交界處、裂紋尖端深度處分別沿與堆芯等高路徑提取溫度。根據(jù)軸向溫度提取線分布情況(見圖8),以Z0處為起點,順時針每隔30°在上述3個位置分別沿一條與堆芯等高的軸線提取路徑上不同時刻的溫度。

圖8 軸向溫度提取線分布情況示意Fig.8 Schematic diagram of distribution of axial temperature extraction lines

繪制溫度分布三維圖如圖9所示(X軸表示環(huán)向角度,(°);Y軸表示軸向高度,換算為無量綱長度:方向為從堆芯下端高度處到堆芯上端高度處;Z軸表示溫度,K)。內(nèi)壁和尖端的冷點在120°和270°(ZN4和Z3處),堆焊層與母材交界處冷點在90°和270°(ZN3和Z3處),所以選定270°(Z3處)為危險位置,Z3處的冷點位置大約在1.9 m高度處。

(a)內(nèi)壁面

(b)堆焊層與母材交界處

(c)裂紋尖端圖9 溫度分布三維圖Fig.9 Three-dimensional diagram of temperature distribution

三進三出RPV模型的冷點位置與二進四出RPV模型的冷點位于DVI兩側(cè)不同[17],三回路模型由于自身進出管口的位置設(shè)計,安注冷點存在一定的偏移。本文已用同樣計算方法驗證過二回路模型,其結(jié)果與文獻[17]中相近,也側(cè)面驗證了模型的有效性。

提取Z3下方1.9 m高度冷點處沿RPV壁面厚度方向溫度變化如圖10、圖11所示,采用無量綱長度作為橫坐標(biāo),0代表內(nèi)壁面,1代表外壁面。

圖10 變化堆芯衰變熱作用下冷點處沿壁厚方向溫度變化Fig.10 Temperature change of the cold spot along the wall thickness direction under the decay heat of the core

圖11 1 MW/m3堆芯衰變熱作用下冷點處沿壁厚方向溫度變化Fig.11 The cold spot temperature changes along the wall thickness direction under the decay heat of the 1 MW/m3 core

在安注水大量注入時(t=106 s),內(nèi)壁溫度比較高,前100 s,在較大的堆芯衰變熱的影響下,堆芯側(cè)壁的流體被加熱,熱量傳遞給吊籃,再傳遞給環(huán)腔內(nèi)的流體,RPV內(nèi)壁一直保持高溫狀態(tài),這導(dǎo)致在安注流量增大時,RPV內(nèi)壁的熱梯度相較1 MW/m3衰變熱時的冷點位置熱梯度大,造成的PTS事故也更嚴重。110 s時,安注流量增大,堆焊層的溫度下降明顯,母材溫度幾乎沒有降低,堆焊層材料導(dǎo)熱系數(shù)約為殼體材料導(dǎo)熱系數(shù)的1/2,堆焊層導(dǎo)熱系數(shù)小,其兩側(cè)的溫度相差較大。110 s之后,安注流速增大,內(nèi)壁溫度下降速度高于外壁,內(nèi)外壁產(chǎn)生一定的溫度梯度,在變形協(xié)調(diào)作用下導(dǎo)致了熱應(yīng)力。但隨著時間逐漸變大,整個壁厚的溫度整體降低。在后續(xù)應(yīng)力分析與斷裂分析中,應(yīng)重點考慮這個危險點(Z3下方1.9 m高度)。

2.3 應(yīng)力分析

沿圖12所示路徑,提取RPV內(nèi)壁等效應(yīng)力的值,并將數(shù)據(jù)繪制成圖(見圖13)?？梢钥闯?t=110 s時,DVI附近兩側(cè)應(yīng)力比較大,在150～200 MPa之間;在t=300 s和500 s時,應(yīng)力有 3個峰值,其中2個峰值較大,1個峰值較小,與溫度三維圖(見圖9)的冷點分布情況相對應(yīng),最大應(yīng)力達374.25 MPa,這與1 MW/m3堆芯衰變熱時最大應(yīng)力322.57 MPa相比,增幅為16.02%;t=1 000 s時峰值應(yīng)力有所降低,但其他位置應(yīng)力有所增大。

圖12 內(nèi)壁應(yīng)力提取路徑Fig.12 Inner wall stress extraction path

圖13 不同時刻環(huán)向應(yīng)力分布Fig.13 Circumferential stress distribution at different times

圖14示出Z3下方1.9 m高度處冷點位置沿壁厚方向應(yīng)力,橫軸作無量綱處理,應(yīng)力提取方向為從內(nèi)壁到外壁,內(nèi)壁為0,外壁為1。

圖14 不同時刻冷點處壁厚方向應(yīng)力分布Fig.14 Stress distribution in the direction of wall thickness at the cold point at different times

由圖14可知,由于內(nèi)壁受低溫安注水冷卻,內(nèi)外壁的溫度梯度大,壁厚方向應(yīng)力呈現(xiàn)內(nèi)拉外壓的情況,且隨著安注時間的延長,外壁溫度也明顯降低,因此外壁應(yīng)力隨時間延長明顯增大,堆焊層與母材交界處存在應(yīng)力不連續(xù)的情況。

2.4 斷裂分析

根據(jù)前述分析可知,Z3下方1.9 m高度處作為整個環(huán)腔區(qū)域內(nèi)的冷點,也是整個環(huán)腔內(nèi)應(yīng)力最大值處,若在此處出現(xiàn)裂紋,則有著最高的斷裂失效風(fēng)險。

文中裂紋深度與長度之比a/c=1/3,假設(shè)裂紋深度a=10 mm,堆焊層內(nèi)深度4 mm,母材內(nèi)深度6 mm。由于假設(shè)裂紋深度較淺,且在小范圍塑性變形內(nèi)應(yīng)力強度因子KI值仍有效[18-19],故用KI作為評估依據(jù),將計算結(jié)果與相應(yīng)條件下的斷裂韌性相比較。斷裂韌性的確定根據(jù)RCC-M附錄ZG 6000中對于核電設(shè)備斷裂評估的相關(guān)規(guī)定,得到斷裂韌性(臨界應(yīng)力強度因子)KIC與裂紋尖端溫度T和材料的韌脆轉(zhuǎn)變參考溫度RTNDT之間的關(guān)系曲線,即KIC=f(T-RTNDT)。其中RTNDT受熱輻照、熱老化和應(yīng)變老化等現(xiàn)象影響。輻照對RTNDT的影響計算式如式(7)所示,其與材料中Cu,P含量有關(guān),也與中子注量f有關(guān)。

ΔRTNDT=[22+556(%Cu-0.08)+2778(%P

-0.008)](f/1019)0.5

(7)

將本文所用的ACP1000模型運行40年各項參數(shù)數(shù)值和母材材料成分參數(shù)數(shù)值代入式(7),算得輻照老化ΔRTNDT=39.667 ℃。被堆焊層覆蓋的母材區(qū)域稱為熱影響區(qū)(HAZ),本模型下,熱老化和應(yīng)變老化引起的ΔRTNDT均為0 ℃,母材韌脆轉(zhuǎn)變溫度主要受輻照影響。經(jīng)計算后,母材斷裂韌性為裂紋尖端溫度的函數(shù),如圖15所示。規(guī)定裂紋尖端應(yīng)力強度因子KI

圖15 母材斷裂韌性曲線Fig.15 Fracture toughness curve of base metal

圖16示出環(huán)腔冷點處深度為10 mm的半橢圓裂紋的尖端應(yīng)力強度因子值。對比1 MW/m3堆芯衰變熱,變化的堆芯衰變熱應(yīng)力強度因子更大,且整體曲線向右偏移。由于106 s后ACC啟動,此時是對內(nèi)壁溫度影響的重要開始時刻,從此刻一直到1 000 s都是主要影響階段,溫度在294～123 ℃(時間500～800 s)之間變化,在此溫度區(qū)間應(yīng)力強度因子較大,兩種模型都出現(xiàn)最大應(yīng)力強度因子,變化堆芯衰變熱相較同時刻1 MW/m3的最大值提高了30.1%。

圖16 應(yīng)力強度因子對比Fig.16 Comparison of stress intensity factors

PTS分析中,按照堆芯衰變熱功率的變化情況來進行計算,可以獲得更加精確的完整性評估結(jié)果,在事故分析中有較高參考價值。

3 結(jié)論

本文在25 cm2小破口失水事故分析中考慮變化的堆芯衰變熱對分析結(jié)果的影響,并與1 MW/m3堆芯衰變熱的結(jié)果相比較,且通過對內(nèi)壁溫度場的提取,確定環(huán)腔內(nèi)壁冷點即溫度最低、應(yīng)力最大的位置,在此處對RPV進行斷裂安全分析。同時,分析了安注水流動軌跡、RPV冷卻路徑、冷卻劑溫度變化等熱工水力問題,得到如下結(jié)論。

(1)在此事故中,初始安注流速小,在DVI下方呈直線下落,后變?yōu)榈妊切?且安注時間越長、安注速度越大,頂角越大。流經(jīng)堆芯的冷卻劑溫度相較1 MW/m3堆芯衰變熱結(jié)果明顯升高,且升溫具有一定的滯后性。

(2)三回路RPV內(nèi)壁安注冷點存在偏移現(xiàn)象,通過提取軸向溫度線,確定環(huán)腔內(nèi)溫度最低點(冷點)位置。與1 MW/m3堆芯衰變熱結(jié)果相比,壁厚方向溫度下降緩慢,變化的堆芯衰變熱對壁面溫度的降低有一定阻礙作用。

(3)沿RPV內(nèi)壁整個環(huán)向的應(yīng)力有3個峰值,與溫度三維圖結(jié)果相對應(yīng),最大應(yīng)力值相較1 MW/m3堆芯衰變熱下的最大應(yīng)力值提高了16.02%。沿冷點厚度方向因內(nèi)部受冷,路徑應(yīng)力呈現(xiàn)內(nèi)拉外壓的情況。

(4)對比1 MW/m3堆芯衰變熱模型,變化的堆芯衰變熱模型的應(yīng)力強度因子明顯提高,且整個曲線向右上方移動,同時刻最大應(yīng)力強度因子提高30.1%,最大應(yīng)力強度因子出現(xiàn)在500～800 s之間。