堆芯衰变热对反应堆压力容器承压热冲击的影响

林 萍,杨森皓,陈旭鹏,银建中,韩志远,谢国山

(1.大连理工大学 化工学院,辽宁大连 116024;2.中国特种设备检测研究院,北京 100029)

0 引言

安全性是核电研究的重中之重,反应堆压力容器(RPV)作为核电站中包容反应堆活性区的重要设备和安全屏障,长期在高温高压高辐照的恶劣工况下服役。当失水事故(LOCA)发生时,反应堆快速停堆,主冷却剂流量降低,应急堆芯冷却系统(ECCS)启动,低温安注水大量注入环腔,产生较大温度梯度,进而产生热应力,与压力的共同作用下产生承压热冲击(PTS)现象[1],使RPV内壁产生裂纹缺陷并持续扩展或使原有缺陷出现扩展[2],甚至有泄漏风险。因而如何更精确地评估PTS现象下RPV的完整性,是目前亟待解决的问题。

PTS工况下,反应堆停堆后由于堆芯内剩余裂变和中子俘获产物衰变产生功率,其功率是一个逐渐下降的过程,开始时下降速度很快,达到一定数值后,下降速度变得缓慢。虽然此值只有稳态功率的百分之几,但仍然对壁面温度的变化有阻碍作用[3-4],直接影响PTS分析结果。马梓淇等[5]在对AP1000型号的RPV在PTS事故下的热应力数值模拟时,考虑了1 MW/m3堆芯热功率的影响。杨森皓等[6]基于三回路RPV模型,进行了PTS瞬态模拟与断裂分析,也设置了1 MW/m3的堆芯衰变热,并对比了传统无堆芯衰变热的计算结果,得到两种PTS事故工况下的衰变热的影响情况。文献[5-6]中均考虑了衰变热的因素,但是所用的衰变热绝对值偏小,且未考虑衰变热是动态变化的过程。

在运行过程中,堆芯区域由于快中子长期辐照导致RPV材料脆化,韧脆转变温度升高[7]。为防止母材的腐蚀和脆化,在堆芯区域RPV内壁有一层厚度3～10 mm的奥氏体不锈钢材料作为堆焊层[8]。堆焊层的相关材料参数与母材材料参数不同,对分析结果有一定影响。

本文在前人研究基础上,基于ACP1000三回路压水堆,对25 cm2小破口失水事故(SBLOCA)工况进行计算;在热工水力计算中,考虑停堆后变化的堆芯衰变热的作用,分析事故下冷却剂流动轨迹、RPV内壁冷却路径,在冷点处假设裂纹进行断裂力学分析;再通过与1 MW/m3堆芯衰变热功率计算结果比较,得到变化堆芯衰变热对同一工况的影响结果。

1 数值计算模型

1.1 物理模型

基于ACP1000三回路压水堆,应用三维绘图软件进行几何建模。外部壳体设计高度为8 000 mm,外径3 200 mm,壁厚200 mm,堆焊层厚度4 mm,壳体内包裹堆芯区域及其他附件,由于其实际结构较复杂,建模过程中进行了相应的简化。如图1所示,其上方交错设置3个进水管和3个出水管,其下方设置2个安注接管(DVI),DVI采用文丘里管形式。低温安注水经由DVI注入RPV环腔内,与腔内原有流体混合,流经下腔室,通过流量分配裙分配流量后,经由堆芯支撑键流入堆芯、带走热量,再从热段出水口流出。将几何模型导入Design Modeler中,进行Slice划分成多个比较规则的形状,有利于后续划分高质量网格,并进行流体填充。进行网格划分时,综合考虑计算精度与计算要求,对DVI等处进行了网格加密处理,最终划分网格数约300万。

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

1.2 堆芯衰变热计算

反应堆停堆后功率包括剩余裂变产生的功率和中子俘获产物产生的功率。停堆后的功率N(t)与停堆前功率N(0)相对变化可表示为式(1),式中右侧第一项是剩余裂变功率,第二项是衰变功率[3]。

(1)

式中,t为停堆时间,s;n(0),n(t)为停堆前和停堆后t时刻的中子密度,中子/cm3;NS(t)为停堆后t时刻的衰变功率,W。

以低浓缩235U作为燃料的轻水堆相对中子密度随时间的变化情况为:

(2)

本文选用的RPV正常运行功率3 180 MW,计算时假设停堆前反应堆已经运行了无限长的时间,由于在计算中忽略了其他俘获产物的衰变功率,故将最终计算结果乘以1.1的安全系数。本文堆芯体积为8.375 64 m3,最终得到2 000 s内单位体积的堆芯衰变热功率如图2所示。

图2 2 000 s内单位体积堆芯衰变热功率变化情况Fig.2 Core decay heat power variation per unit volume in 2 000 s

1.3 控制方程

环腔内流体假设为粘性不可压缩湍流模型,其基本方程为Reynolds时均方程,选择Realizablek-ε两方程模型[9],此模型Cμ是与ε有关的公式。耗散率传输方程相较标准模型在生成项处不包含Pk,减少项中不具有奇异性。Realizablek-ε湍流模型控制方程[10]如下:

(3)

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

流体域计算的温度、压力等参数与壁面的耦合边界上,要满足温度和热流密度的连续性条件:

(4)

式中,Tfp,Tsp为交界面两侧流体、固体的网格节点温度;qf,qs为流体、固体两侧的热流密度。

冷却剂与RPV接触,热量传递主要形式为对流换热及热传导,对流换热边界条件为:

(5)

式中,h为流体和固体间的表面对流换热系数,W/(m2·℃);TS,Tf为RPV壁面温度、流体温度,℃。

固体传热方程为:

(6)

式中,ρ为材料密度,kg/m3;Cp为材料的定压比热容,J/(kg·℃);λ为材料的导热系数,W/(m·℃);qv为堆芯单位体积释热率,与堆芯热功率有关,W/m3。

1.4 求解设置

1.4.1 事故工况

本文计算的PTS工况是25 cm2小破口失水事故工况,该工况下的输入温度、压力条件见图3,在冷段进水口设置流量、温度入口,热段作为压力出口。DVI接口安注水温度为10 ℃,小破口失水事故中ECCS以压力作为触发条件,t=6.5 s触发堆芯补水箱(CMT),在t=106 s时内压下降到蓄压安注箱(ACC)的触发条件,具体参数见表1。

图3 25 cm2小破口失水事故输入条件Fig.3 Input conditions for 25 cm2 small break loss of coolant accident

表1 冷却剂注入各项参数Tab.1 Coolant injection parameters

1.4.2 条件设置

基于CFX平台进行流体仿真计算,堆芯简化为多孔模型[11-12],将冷却剂、堆芯、RPV外壳、堆焊层划分为4个区域,通过设置Fluid-Porous型、Fluid-Solid型和Solid-Solid型的交界面,实现不同区域之间的耦合。

RPV外壳和堆焊层两者均设置为固体域,RPV壳体、接管和主要结构材料为16MnNiMo(16MND5),堆焊层材料为奥氏体不锈钢[13],材料物性参数分别如表2、表3所示。RPV外壁绝热,材料计算模型应用双线性等向强化的弹塑性模型。冷却剂为流体域,介质为水,水的物性采用IAPWS-IF97标准计算,初始采用温度300 ℃、压力15.267 MPa时的值,在计算中考虑浮力的影响,水的物性随温度、压力变化。堆芯材料为二氧化铀和锆锡合金陶瓷芯块,堆芯为多孔域,流动损失系数0.39 m-1,横向损失系数10 m-1,设置的热源项为随时间变化的函数(见图2)。

表2 16MND5材料参数Tab.2 Material parameters of 16MND5

表3 堆焊层材料参数Tab.3 Material parameters of surfacing layer

2 结果分析

2.1 流动混合分析

由于安注水在环腔内与原有冷却剂混合的同时伴随着热量的交换,因此,通过提取环腔内冷却剂温度变化情况来分析安注水注入环腔后的流动轨迹。安注水注入环腔后的温度变化情况如图4所示。

图4 安注水在环腔内流动混合情况Fig.4 The mixing condition of water flow in the annular cavity

当只有CMT启动时,安注水流速较小,安注外轮廓轨迹近似为一条直线;随着安注时间变长,安注流体与环腔内冷却剂在环腔下部混合,且逐渐向整个环腔的周向区域扩散。t>106 s时,ACC启动,安注流速增大,安注水流动轨迹近似为等腰三角形;随着安注时间的增加,等腰三角形的顶角更大,环腔内安注水大范围扩散,环腔内混合水温开始逐渐降低,二者混合逐渐均匀。

安注水流经环腔与冷却剂混合后,经由下腔室流入堆芯。从堆芯中部横截面位置的温度分布图(见图5(a))可以看出,堆芯内水温最高的区域集中在中心轴线附近,所以在堆芯中心线处沿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时,伴随着前期高热功率的累积,堆芯温度仍然在升高,但堆芯周围冷却剂逐渐循环带走下方部分热量,A,B点温差最大为190.77 K左右,与理论流体吸收热量计算公式Q=cmΔT所计算出的温差相比误差为6.41%,误差在允许范围内[14]。时间递进到110 s,ACC已经启动,大量低温安注水注入RPV,堆芯冷却速率加快,堆芯温度较之前时刻降低,安注水流量增加,冷却剂温度降低,到达B点可换热的冷却液多,B点温度下降。同条件下,计算1 MW/m3堆芯衰变热作用下,B点相较A点温度升高了30 K左右。

冷却剂在流动过程中受压力、浮力以及堆芯多孔域的影响,经过堆芯的水流动轨迹存在一定的偏移现象,各处流体被流经堆芯带走热量的路径不可确定,所以不是越接近B点堆芯温度越高。

安注冷却剂在RPV环腔内的流动情况会影响RPV内壁的温度分布,进而在内壁产生较大热应力,所以在进行反应堆压力容器的完整性评估时,要根据安注冷却剂流动轨迹,分析RPV内壁上温度的变化。

2.2 传热分析

ECCS启动后,RPV内壁面温度场分布如图7所示。低速安注初期,DVI下方有一条竖直的羽流区域[15-16],安注时间增大,羽流区域发生偏移,偏向DVI两侧。安注流量增大后,羽流区域更偏向于两侧,此现象导致RPV环腔周向产生较大的温度梯度,对环焊缝、缺陷强度和完整性评估有重要影响。

图7 RPV内壁温度变化Fig.7 Inner wall temperature variation of RPV

进行完整性评估时,确定安注冷却过程冷点的位置至关重要,冷点处温度梯度大,产生的热应力也大,若此处存在缺陷,则最容易发生裂纹扩展和失效现象。为确定冷点位置,在RPV内壁、堆焊层与母材交界处、裂纹尖端深度处分别沿与堆芯等高路径提取温度。根据轴向温度提取线分布情况(见图8),以Z0处为起点,顺时针每隔30°在上述3个位置分别沿一条与堆芯等高的轴线提取路径上不同时刻的温度。

图8 轴向温度提取线分布情况示意Fig.8 Schematic diagram of distribution of axial temperature extraction lines

绘制温度分布三维图如图9所示(X轴表示环向角度,(°);Y轴表示轴向高度,换算为无量纲长度:方向为从堆芯下端高度处到堆芯上端高度处;Z轴表示温度,K)。内壁和尖端的冷点在120°和270°(ZN4和Z3处),堆焊层与母材交界处冷点在90°和270°(ZN3和Z3处),所以选定270°(Z3处)为危险位置,Z3处的冷点位置大约在1.9 m高度处。

(a)内壁面

(b)堆焊层与母材交界处

(c)裂纹尖端图9 温度分布三维图Fig.9 Three-dimensional diagram of temperature distribution

三进三出RPV模型的冷点位置与二进四出RPV模型的冷点位于DVI两侧不同[17],三回路模型由于自身进出管口的位置设计,安注冷点存在一定的偏移。本文已用同样计算方法验证过二回路模型,其结果与文献[17]中相近,也侧面验证了模型的有效性。

提取Z3下方1.9 m高度冷点处沿RPV壁面厚度方向温度变化如图10、图11所示,采用无量纲长度作为横坐标,0代表内壁面,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),内壁温度比较高,前100 s,在较大的堆芯衰变热的影响下,堆芯侧壁的流体被加热,热量传递给吊篮,再传递给环腔内的流体,RPV内壁一直保持高温状态,这导致在安注流量增大时,RPV内壁的热梯度相较1 MW/m3衰变热时的冷点位置热梯度大,造成的PTS事故也更严重。110 s时,安注流量增大,堆焊层的温度下降明显,母材温度几乎没有降低,堆焊层材料导热系数约为壳体材料导热系数的1/2,堆焊层导热系数小,其两侧的温度相差较大。110 s之后,安注流速增大,内壁温度下降速度高于外壁,内外壁产生一定的温度梯度,在变形协调作用下导致了热应力。但随着时间逐渐变大,整个壁厚的温度整体降低。在后续应力分析与断裂分析中,应重点考虑这个危险点(Z3下方1.9 m高度)。

2.3 应力分析

沿图12所示路径,提取RPV内壁等效应力的值,并将数据绘制成图(见图13)。可以看出,t=110 s时,DVI附近两侧应力比较大,在150～200 MPa之间;在t=300 s和500 s时,应力有 3个峰值,其中2个峰值较大,1个峰值较小,与温度三维图(见图9)的冷点分布情况相对应,最大应力达374.25 MPa,这与1 MW/m3堆芯衰变热时最大应力322.57 MPa相比,增幅为16.02%;t=1 000 s时峰值应力有所降低,但其他位置应力有所增大。

图12 内壁应力提取路径Fig.12 Inner wall stress extraction path

图13 不同时刻环向应力分布Fig.13 Circumferential stress distribution at different times

图14示出Z3下方1.9 m高度处冷点位置沿壁厚方向应力,横轴作无量纲处理,应力提取方向为从内壁到外壁,内壁为0,外壁为1。

图14 不同时刻冷点处壁厚方向应力分布Fig.14 Stress distribution in the direction of wall thickness at the cold point at different times

由图14可知,由于内壁受低温安注水冷却,内外壁的温度梯度大,壁厚方向应力呈现内拉外压的情况,且随着安注时间的延长,外壁温度也明显降低,因此外壁应力随时间延长明显增大,堆焊层与母材交界处存在应力不连续的情况。

2.4 断裂分析

根据前述分析可知,Z3下方1.9 m高度处作为整个环腔区域内的冷点,也是整个环腔内应力最大值处,若在此处出现裂纹,则有着最高的断裂失效风险。

文中裂纹深度与长度之比a/c=1/3,假设裂纹深度a=10 mm,堆焊层内深度4 mm,母材内深度6 mm。由于假设裂纹深度较浅,且在小范围塑性变形内应力强度因子KI值仍有效[18-19],故用KI作为评估依据,将计算结果与相应条件下的断裂韧性相比较。断裂韧性的确定根据RCC-M附录ZG 6000中对于核电设备断裂评估的相关规定,得到断裂韧性(临界应力强度因子)KIC与裂纹尖端温度T和材料的韧脆转变参考温度RTNDT之间的关系曲线,即KIC=f(T-RTNDT)。其中RTNDT受热辐照、热老化和应变老化等现象影响。辐照对RTNDT的影响计算式如式(7)所示,其与材料中Cu,P含量有关,也与中子注量f有关。

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

-0.008)](f/1019)0.5

(7)

将本文所用的ACP1000模型运行40年各项参数数值和母材材料成分参数数值代入式(7),算得辐照老化ΔRTNDT=39.667 ℃。被堆焊层覆盖的母材区域称为热影响区(HAZ),本模型下,热老化和应变老化引起的ΔRTNDT均为0 ℃,母材韧脆转变温度主要受辐照影响。经计算后,母材断裂韧性为裂纹尖端温度的函数,如图15所示。规定裂纹尖端应力强度因子KI

图15 母材断裂韧性曲线Fig.15 Fracture toughness curve of base metal

图16示出环腔冷点处深度为10 mm的半椭圆裂纹的尖端应力强度因子值。对比1 MW/m3堆芯衰变热,变化的堆芯衰变热应力强度因子更大,且整体曲线向右偏移。由于106 s后ACC启动,此时是对内壁温度影响的重要开始时刻,从此刻一直到1 000 s都是主要影响阶段,温度在294～123 ℃(时间500～800 s)之间变化,在此温度区间应力强度因子较大,两种模型都出现最大应力强度因子,变化堆芯衰变热相较同时刻1 MW/m3的最大值提高了30.1%。

图16 应力强度因子对比Fig.16 Comparison of stress intensity factors

PTS分析中,按照堆芯衰变热功率的变化情况来进行计算,可以获得更加精确的完整性评估结果,在事故分析中有较高参考价值。

3 结论

本文在25 cm2小破口失水事故分析中考虑变化的堆芯衰变热对分析结果的影响,并与1 MW/m3堆芯衰变热的结果相比较,且通过对内壁温度场的提取,确定环腔内壁冷点即温度最低、应力最大的位置,在此处对RPV进行断裂安全分析。同时,分析了安注水流动轨迹、RPV冷却路径、冷却剂温度变化等热工水力问题,得到如下结论。

(1)在此事故中,初始安注流速小,在DVI下方呈直线下落,后变为等腰三角形,且安注时间越长、安注速度越大,顶角越大。流经堆芯的冷却剂温度相较1 MW/m3堆芯衰变热结果明显升高,且升温具有一定的滞后性。

(2)三回路RPV内壁安注冷点存在偏移现象,通过提取轴向温度线,确定环腔内温度最低点(冷点)位置。与1 MW/m3堆芯衰变热结果相比,壁厚方向温度下降缓慢,变化的堆芯衰变热对壁面温度的降低有一定阻碍作用。

(3)沿RPV内壁整个环向的应力有3个峰值,与温度三维图结果相对应,最大应力值相较1 MW/m3堆芯衰变热下的最大应力值提高了16.02%。沿冷点厚度方向因内部受冷,路径应力呈现内拉外压的情况。

(4)对比1 MW/m3堆芯衰变热模型,变化的堆芯衰变热模型的应力强度因子明显提高,且整个曲线向右上方移动,同时刻最大应力强度因子提高30.1%,最大应力强度因子出现在500～800 s之间。