陳明亞，耿昌金，王威強(qiáng)，高紅波，彭群家，師金華

(1.蘇州熱工研究院有限公司，江蘇蘇州 215004;2.山東大學(xué) 機(jī)械工程學(xué)院，濟(jì)南 250061)

0 引言

反應(yīng)堆壓力容器(RPV)是核安全一級(jí)部件，在服役過(guò)程中，由于受到中子輻照的影響，材料性能將會(huì)逐漸劣化，具體表現(xiàn)為強(qiáng)度增加、塑韌性下降[1-2]。同時(shí)，RPV在制造、安裝、服役過(guò)程中又會(huì)出現(xiàn)裂紋類缺陷，這些因素都將對(duì)核電站的安全運(yùn)行產(chǎn)生嚴(yán)重的影響[3]。因此，在設(shè)計(jì)階段需要進(jìn)行含假想裂紋的斷裂力學(xué)安全性能評(píng)估，并且在運(yùn)行過(guò)程中若發(fā)生超出設(shè)計(jì)運(yùn)行壓力-溫度限值曲線(P-T曲線)時(shí)，也需要進(jìn)行含假想裂紋的斷裂力學(xué)安全性能評(píng)估[4-5]。

在RPV輻照脆化工程評(píng)估中，需要考慮幾何信息(含參數(shù)敏感性分析)、瞬態(tài)參數(shù)(含大量設(shè)計(jì)瞬態(tài)及實(shí)際運(yùn)行瞬態(tài))、材料(含輻照效應(yīng)引起的變化)、斷裂評(píng)估準(zhǔn)則等影響參數(shù)[6-7]。同時(shí)，斷裂參量(應(yīng)力強(qiáng)度因子等)計(jì)算過(guò)程復(fù)雜，難以滿足工程上快速評(píng)價(jià)的應(yīng)用需求[8-10]。為此，進(jìn)行RPV輻照脆化力學(xué)評(píng)估專用分析平臺(tái)的開(kāi)發(fā)十分必要。有限元軟件Ansys的APDL編程語(yǔ)言(Ansys Parametric Design Language，即Ansys參數(shù)化設(shè)計(jì)語(yǔ)言)提供了廣泛的二次開(kāi)發(fā)潛能，APDL語(yǔ)言是一種類似Fortran的解釋性語(yǔ)言，提供一般程序語(yǔ)言的功能，如參數(shù)、宏、標(biāo)量、向量及矩陣運(yùn)算、分支、循環(huán)、重復(fù)以及訪問(wèn)Ansys有限元數(shù)據(jù)庫(kù)等，另外還提供簡(jiǎn)單界面定制功能，實(shí)現(xiàn)參數(shù)交互輸入、消息機(jī)制、界面驅(qū)動(dòng)和運(yùn)行應(yīng)用程序等?？苫贏PDL語(yǔ)言開(kāi)發(fā)RPV輻照脆化評(píng)估參數(shù)化專用模塊，參數(shù)化專用模塊可集成評(píng)估模型基本信息輸入、溫度場(chǎng)計(jì)算、應(yīng)力場(chǎng)計(jì)算、斷裂參量計(jì)算、依據(jù)RCC-M規(guī)范進(jìn)行安全評(píng)估等方面的分析能力。開(kāi)發(fā)參數(shù)化專用模塊可規(guī)范計(jì)算過(guò)程，避免了人因干擾，滿足工程上的快速、準(zhǔn)確的安全評(píng)估要求。

本文首先梳理RCC-M規(guī)范對(duì)RPV輻照脆化評(píng)估的分析要求，然后基于大型有限元軟件Ansys平臺(tái)開(kāi)發(fā)RPV輻照脆化評(píng)估參數(shù)化專用模塊，最后通過(guò)與某核電廠原設(shè)計(jì)報(bào)告中分析結(jié)果進(jìn)行對(duì)比，驗(yàn)證參數(shù)化專用模塊的可靠性。

1 分析規(guī)范要求

1.1 分析模型與評(píng)估準(zhǔn)則

RPV堆芯段受到的快中子注量最大，某RPV堆芯筒體段的幾何結(jié)構(gòu)如圖1所示。分析案例中，筒體內(nèi)徑Ri為1 994.5 mm、筒體母材壁厚t為200 mm、筒體內(nèi)表面堆焊層厚度tcladding為7.5 mm[11-12]。

圖1 RPV堆芯段分析模型Fig.1 RPV core segment model

參考RCC-M規(guī)范確定不同工況下的評(píng)估準(zhǔn)則，其中正常與擾動(dòng)工況下安全性能評(píng)價(jià)準(zhǔn)則如式(1)所示。應(yīng)力強(qiáng)度因子KⅠ是根據(jù)所選擇的分析時(shí)刻有關(guān)的薄膜應(yīng)力強(qiáng)度因子和彎曲應(yīng)力強(qiáng)度因子之和，一次薄膜應(yīng)力和一次彎曲應(yīng)力引起的應(yīng)力強(qiáng)度因子必須乘以評(píng)定系數(shù)2。要求輸出整個(gè)瞬態(tài)中KⅠ/KⅠR的最大值及相應(yīng)的時(shí)刻t，RCC-M規(guī)范要求KⅠ/KⅠR的最大值小于1。

KⅠ=2KⅠm+KⅠt≤KⅠR

(1)

式中,KⅠm為總體一次應(yīng)力強(qiáng)度因子,MPa·m1/2；KⅠt為由溫差產(chǎn)生的應(yīng)力強(qiáng)度因子， MPa·m1/2；KⅠR為材料參考斷裂韌度，MPa·m1/2。

本文對(duì)緊急、事故和水壓試驗(yàn)的評(píng)價(jià)準(zhǔn)則不再贅述，后續(xù)內(nèi)容也僅針對(duì)正常與擾動(dòng)工況的應(yīng)用展開(kāi)論述。

1.2 材料斷裂韌性分析

堆芯材料基本物理性能參考RCC-M規(guī)范中的數(shù)據(jù)，材料參考斷裂韌性KⅠC和KⅠa分別如式(2)(3)所示，取KⅠC和KⅠa的最小值。

KⅠC=min{36.5+3.1exp[0.36(T-RTNDT

+55.5)],220}

(2)

KⅠa=min{29.43+1.355exp[0.0261(T-RTNDT

+88.9)],220}

(3)

式中,T為溫度參數(shù),℃;RTNDT為材料的韌脆轉(zhuǎn)變溫度,℃。

RPV材料韌脆轉(zhuǎn)變溫度RTNDT計(jì)算方法如下:

RTNDT=RTNDT(i)+ΔRTNDT

(4)

式中，RTNDT(i)為初始韌脆轉(zhuǎn)變溫度，℃；ΔRTNDT為韌脆轉(zhuǎn)變溫度增量，℃。

采用RCC-M規(guī)范推薦的預(yù)測(cè)公式，計(jì)算ΔRTNDT：

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

×(w(P)-0.008)](f/1019)1/2

(5)

式中，f為中子注量，n/cm2;w(Cu)，w(P)分別為RPV材料中銅元素和磷元素的質(zhì)量百分含量，具體應(yīng)用中，如果w(Cu)<0.08，則w(Cu)取0.08，如果w(P)<0.008，則w(P)取0.008。

1.3 應(yīng)力強(qiáng)度因子計(jì)算

通過(guò)所考慮的缺陷尺寸和分析工況下的應(yīng)力來(lái)確定應(yīng)力強(qiáng)度因子KⅠ，具體流程如下。

(1)確定所分析的工況及該工況的瞬時(shí)應(yīng)力分布(考慮所有的施加載荷)。分析時(shí)考慮沿著所確定分析截面的支承部位，在假定缺陷平面上的正應(yīng)力分布。用離內(nèi)壁面的距離x表示支承部位上的每一點(diǎn)，正應(yīng)力(σ(x))的分布由變量x除以截面厚度t的多項(xiàng)式表示:

σ(x)=σ0+σ1(x/t)+σ2(x/t)2+σ3(x/t)3

+σ4(x/t)4

(6)

(2)結(jié)合影響函數(shù)確定KⅠ，影響函數(shù)的系數(shù)分別由i0，i1，i2,i3表示，是缺陷(裂紋)幾何形狀、假定裂紋所在區(qū)域及比值a/t的函數(shù)(a為裂紋深度)。式(7)給出了結(jié)合影響函數(shù)計(jì)算KⅠ的表達(dá)式。

KⅠ=(πa)1/2[σ0i0+σ1(a/t)i1+σ2(a/t)2i2

+σ3(a/t)3i3]

(7)

表1列出了裂紋形狀系數(shù)a/b=1/3的半橢圓裂紋的影響函數(shù)值i0，i1,i2,i3，表中的影響函數(shù)值適用于R/t≤10的空心圓柱體模型(對(duì)應(yīng)于本文中的RPV結(jié)構(gòu))。

表1 半橢圓裂紋(a/b=1/3)影響函的系數(shù)取值Tab.1 Values of influence function for semi-elliptical crack(a/b=1/3)

2 軟件參數(shù)化評(píng)估

2.1 軟件編程框架

參數(shù)化專用模塊的基本編程框架如圖2所示，參數(shù)化專用模塊界面如圖3所示，專用模塊新增按鈕功能如表2所示。此程序設(shè)有RPV信息輸入和瞬態(tài)信息輸入兩個(gè)輸入模塊，包含溫度場(chǎng)計(jì)算、應(yīng)力場(chǎng)計(jì)算和結(jié)構(gòu)安全性能評(píng)估3個(gè)分析模塊。

表2 基本按鈕功能介紹Tab.2 Introduction of basic functions of buttons

圖2 參數(shù)化專用模塊基本編程框架Fig.2 Basic programming framework of the special plug-in module

圖3 參數(shù)化專用模塊界面Fig.3 Interface of the special plug-in module

2.2 基本輸入信息

2.2.1 堆芯模型及假想缺陷信息

通過(guò)“RPV_INFOR”按鈕設(shè)計(jì)堆芯模型相關(guān)信息。軟件提供了堆芯模型的基本信息輸入界面，如圖4所示。

R為堆芯筒體內(nèi)徑，m；Tbase為堆芯筒體母材厚度，m；Tclad為堆芯筒體堆焊層厚度，m；RTNDTi為堆芯筒體母材初始韌脆轉(zhuǎn)變溫度,℃；WCu為堆芯筒體母材中Cu元素的質(zhì)量百分含量(%)；WP為堆芯筒體母材中P元素的質(zhì)量百分含量(%)；F為堆芯筒體的快中子注量(n/cm2，能量大于1 MeV的快中子)圖4 堆芯模型基本信息輸入界面Fig.4 Basic information input interface of the core segment

所考慮的基準(zhǔn)缺陷為一平面型半橢圓表面裂紋，該裂紋位于應(yīng)力最大的表面上，并假設(shè)此缺陷平面垂直于最大主應(yīng)力方向，如圖1所示。在正常與擾動(dòng)工況下，基準(zhǔn)缺陷的深度尺寸a為1/4容器厚度t。分析中所用的堆芯筒體段有限元模型如圖5所示，共包含1 920個(gè)有限元單元。

圖5 堆芯有限元模型Fig.5 Finite element model of the core segment

2.2.2 分析瞬態(tài)信息

通過(guò)“TRANSIENT_PARA”按鈕輸入瞬態(tài)信息，輸入界面如圖6所示。

Ini_T 為瞬態(tài)的初始溫度，℃；Nstep為瞬態(tài)的載荷步，℃；TR為瞬態(tài)的總時(shí)間，s圖6 瞬態(tài)相關(guān)信息輸入界面Fig.6 Input interface for the transient related information

瞬態(tài)具體的壓力、溫度隨時(shí)間的變化關(guān)系由原設(shè)計(jì)報(bào)告提取，或由運(yùn)行監(jiān)測(cè)報(bào)告中提取，并按照指定的格式編寫.inp或.txt格式的數(shù)據(jù)文件。

2.3 有限元數(shù)值仿真與安全評(píng)價(jià)

2.3.1 有限元數(shù)值仿真

通過(guò)“TEMPERATURE_CAL”和“STRESS_CAL”兩個(gè)界面按鈕調(diào)用Ansys本身模塊進(jìn)行溫度場(chǎng)與應(yīng)力場(chǎng)計(jì)算。本軟件使用的是Plan 77，該單元為二維八節(jié)點(diǎn)四邊形單元。

在正常與擾動(dòng)工況，RPV內(nèi)表面與流體之間換熱系數(shù)保守的取為無(wú)窮大，外表面取為絕熱邊界條件。有限元數(shù)值分析中采用二維軸對(duì)稱模型，分別耦合圖1中上下截面的軸向位移。

2.3.2 Ansys矩陣運(yùn)算

矩陣運(yùn)算是一種數(shù)組參數(shù)之間的數(shù)學(xué)運(yùn)算，例如矩陣乘法、計(jì)算轉(zhuǎn)置矩陣、求解聯(lián)立方程組等。 Ansys軟件對(duì)兩個(gè)輸入數(shù)組參數(shù)矩陣進(jìn)行矩陣運(yùn)算，輸出一個(gè)數(shù)組參數(shù)矩陣。矩陣運(yùn)算包括：(1)矩陣相乘；(2)求解聯(lián)立方程組；(3)對(duì)矩陣中的某個(gè)指定向量排序；(4)計(jì)算兩個(gè)向量之間的協(xié)方差；(5)計(jì)算兩個(gè)向量之間的相關(guān)性等。

本軟件使用的矩陣運(yùn)算的主要命令包括：*MOPER或Utility Menu>Parameters>Array Operations>Matrix Operations。

*MOPER命令可以求解聯(lián)立方程組系數(shù)的矩陣，根據(jù)評(píng)估點(diǎn)的坐標(biāo)，將式(5)調(diào)整成方程組的形式，如式(8)所示。

ai1X1+ai2X2+…+ainXn=bi

(8)

式(8)中位置矩陣X，應(yīng)力結(jié)果矩陣b，通過(guò)Ansys軟件可以求解出系數(shù)矩陣A，命令如下所示：

*MOPER,b,X,SOLV,A

2.3.3 安全性能評(píng)價(jià)

通過(guò)“SAFETY_ASSESS”按鈕，可自動(dòng)進(jìn)行依據(jù)RCC-M規(guī)范的安全性能評(píng)價(jià)。

3 軟件展示與驗(yàn)證說(shuō)明

通過(guò)某一瞬態(tài)的安全性能評(píng)價(jià)展示參數(shù)化專用模塊的使用過(guò)程，其主要運(yùn)用過(guò)程說(shuō)明如下。

(1)Step 1：在Ansys軟件的“Command Prompt”輸入框中輸入“SNPI_Button”運(yùn)行相關(guān)的宏文件。

(2)Step 2：運(yùn)行“RPV_INFOR”按鈕，輸入RPV的基本輸入信息，堆芯區(qū)域母材為16MND5(法國(guó)RCC-M規(guī)范中的M2111材料)，堆焊層材料為E309L和E308L。

(3)Step 3：運(yùn)行“TRANSIENT_PARA”按鈕，定義瞬態(tài)的相關(guān)信息，分析案例的瞬態(tài)信息如圖7所示，根據(jù)瞬態(tài)特性編寫瞬態(tài)輸入信息文件“Tr_Tem_Pre.inp”(瞬態(tài)具體的壓力、溫度隨時(shí)間的變化關(guān)系按照指定的格式編寫.inp文件，如圖8所示)。

圖7 分析案例的瞬態(tài)信息Fig.7 Transient information of the studied case

圖8 壓力、溫度隨時(shí)間變化關(guān)系的.inp輸入文件Fig.8 Time-depended data of pressure and temperature in the .inp file

(4)Step 4：運(yùn)行“TEMPERATURE_CAL”按鈕，進(jìn)行瞬態(tài)溫度場(chǎng)計(jì)算。

(5)Step 5：運(yùn)行“STRESS_CAL”按鈕，進(jìn)行瞬態(tài)應(yīng)力場(chǎng)計(jì)算。

(6)Step 6：運(yùn)行“SAFETY_ASSESS”按鈕，進(jìn)行安全性能評(píng)價(jià)。RPV堆芯結(jié)構(gòu)內(nèi)表面缺陷分析數(shù)據(jù)輸出數(shù)據(jù)示例如圖9所示，安全分析輸出數(shù)據(jù)如表3所示(Office Excel的數(shù)據(jù)格式)。

表3 軸向內(nèi)表面缺陷安全分析輸出數(shù)據(jù)Tab.3 Output data of safety assessment for internal axial surface crack

圖9 軸向內(nèi)表面缺陷分析數(shù)據(jù)輸出Fig.9 Data output for internal axial surface crack

專用模塊可以獲得瞬態(tài)中結(jié)構(gòu)最危險(xiǎn)的時(shí)刻及結(jié)構(gòu)的安全性能。該瞬態(tài)中KⅠ/KⅠR的最大值為0.78，滿足RCC-M規(guī)范要求，獲得最大值對(duì)應(yīng)的瞬態(tài)時(shí)間為52 000 s。

軟件可靠性驗(yàn)證結(jié)果表明，參數(shù)化專用模塊的分析結(jié)果與某核電廠原設(shè)計(jì)報(bào)告中瞬態(tài)的分析結(jié)果偏差可控制在3%左右，依據(jù)IAEA的工程評(píng)估對(duì)比原則[13-14]，認(rèn)為本軟件可獲取與原設(shè)計(jì)報(bào)告一致的分析結(jié)果。

4 結(jié)語(yǔ)

基于Ansys軟件自身的APDL語(yǔ)言開(kāi)發(fā)了RPV輻照脆化評(píng)估參數(shù)化專用模塊，并以某正常與擾動(dòng)工況為例，進(jìn)行了開(kāi)發(fā)的專用模塊應(yīng)用展示。開(kāi)發(fā)的專用軟件模塊規(guī)范了計(jì)算過(guò)程，避免了人因干擾，可滿足工程上快速、準(zhǔn)確的安全評(píng)估要求，并且對(duì)比驗(yàn)證表明，參數(shù)化專用模塊的分析結(jié)果與某核電廠的原設(shè)計(jì)報(bào)告分析結(jié)果一致。