王大勝，周 強，毛劍峰

(1.深圳中廣核工程設計有限公司，廣東深圳 518124；2.浙江工業(yè)大學 化工機械設計研究所，杭州 310023)

0 引言

在堆芯熔毀的嚴重事故(SA)中，先進的核電站(NPP)廣泛采用 “容器內滯留(IVR)”策略來維持反應堆壓力容器(RPV)的安全[1]。事實上，自1996年以來，美國核管理委員會(NRC)已將IVR措施認證為SA管理的標準措施。在大多數(shù)情況下，IVR措施是在RPV的外壁上提供長期的水冷卻，因此在SA期間無需任何人為行動，就可以去除衰變熱。IVR措施實施過程中，外部水冷卻是事故管理的最重要特征，如圖1[2]所示。在IVR的傳統(tǒng)概念中，RPV的下封頭(LH)在堆芯碎片到達內部之前完全淹沒在水浸中[3]。因此，在RPV內形成的熔池溫度約為1 327 ℃，而容器外壁的溫度接近150 ℃[4]。根據(jù)事故期間RPV臨界熱流(CHF)的基本假設，可以從先前的研究中得知，RPV的安全性可以在規(guī)定的時間內得到保證[5-7]。在2011年日本福島核事故之前，上述常規(guī)IVR概念并未受到嚴重挑戰(zhàn)。但是，當冷卻水突然注入后，RPV結構內外溫差急劇增大，加上熱沖擊作用，RPV外表面產(chǎn)生了多處裂紋，如圖2，3[8]所示，事故后果表明，按照無裂紋RPV結構的安全評估結果受到了嚴重的挑戰(zhàn)[9]。實際上，標準事故設計中沒有考慮帶裂紋的RPV結構，因而類似福島核事故完全超出了傳統(tǒng)事故條件下的設計和安全評估范圍。因此，在嚴重事故條件的72 h內，有裂紋的RPV能否保持壓力邊界的完整性是防止核泄漏的關鍵問題之一。

圖2 失效后的反應堆壓力容器Fig.2 The surface of failed RPV

歐盟核電堆芯熔毀試驗項目(EC-FOREVER)測試提供了豐富的數(shù)據(jù)源，部分驗證了RPV的斷裂模型和歐盟核電行業(yè)結構安全評估規(guī)范，其RPV破裂如圖3所示。EC-FOREVER的項目總協(xié)調在德國羅森多夫進行，以WILLSCHüTZ等[10]都對每項測試進行背靠背地模擬計算及分析，主要采用Ansys多物理耦合方法，建立2D的RPV結構對主要IVR測試項目進行了熱-流-固耦合分析。國內MAO等[11]研究者利用非線性多物理場有限元法,計算分析了IVR條件下RPV的蠕變失效。KOUNDY等[12]指出，測量和計算的最大位移值都出現(xiàn)在失效破裂處，Ansys模擬預測的RPV破裂點與試驗測量值非常一致。對于IVR下RPV結構分析，一般需要同時考慮熱模塊和結構模塊的耦合計算，才能較好地預測RPV 失效模式、位置和時間。KULKARNI等[13]利用38 kW的恒定功率加熱器對RPV進行堆芯熔融試驗，容器內充入25 bar的氬氣，發(fā)現(xiàn)RPV材料的初始熔化溫度略小于焊縫熔化溫度，利用數(shù)值方法分析了2D軸對稱RPV失效破裂，發(fā)現(xiàn)裂紋水平沿壁厚法向分布，并分析了單個裂紋行為[14]。ABENDROTH等[15]計算分析了包括1/1(真實比例)和1/10(縮比模型)的3D裂紋尺寸及其位置變化對RPV破裂行為影響。本文給出RPV在嚴重事故條件下的斷裂力學計算結果，包括與溫度相關的應力(應變)分布。由于外壁溫度低于蠕變溫度范圍，沒有考慮蠕變影響，但是建立了含裂紋下封頭的有限元模型，用于分析研究J積分、溫度及應力變化規(guī)律。

圖3 裂紋和蠕變孔隙破壞示意圖Fig.3 The schematic diagram of creep porosity and crack

1 理論背景與建模

1.1 熱應力體J積分的理論模型

為了簡單起見，先給出適用于超彈性材料的下列公式。

(1)

(2)

式中，σij為熱-機引起的勢函數(shù),表示應力張量；W為應變能密度；εij為應變張量。

實際上，上述公式也可用于塑性增量理論。在該理論中，禁止局部卸載，應力空間中的所有加載路徑應保持徑向，以便主應力的比率保持不變。塑性區(qū)應封閉在積分路徑(或積分區(qū)域)內，并假設裂紋面沒有載荷，也不存在體載荷。

為了說明3DJ積分，繪制了具有彎曲裂紋前緣的3D平面裂紋,如圖4所示，圖中定義了一個局部裂紋尖端坐標系，其中ξ1軸垂直于裂紋擴展方向上的裂紋前緣，而ξ2軸垂直于裂紋面，ξ3軸與裂紋前緣相切，s是局部裂紋前緣坐標。實際上，在計算3DJ積分時，通常使用兩種典型的J積分公式，第一種是路徑積分和面積積分的組合，即：

圖4 3D積分路徑/面積Fig.4 The 3D integral path/area

(3)

關于虛擬裂紋擴展(VCE)的另一個公式[16]，采用qk的3D坐標函數(shù)和體積積分域來求解(見圖5)：

(4)

在上述兩個方程中，指數(shù)i，j，k在坐標ξ1,ξ2,ξ3上變化。式(4)中，q1是ξ3和裂紋尖端徑向距離(r，q2=q3=0)的函數(shù)，q1可表示為:

(5)

具體講，積分域是具有長度δs和半徑r0的圓柱形體積，如圖5所示。由此產(chǎn)生的虛擬裂紋前緣擴展(VCFE)可描述為一個平面拋物線形區(qū)域，其大小可表示為：

圖5 具有彎曲裂紋前緣的3D積分體積Fig.5 The 3D integral volume of curved crack front edge

(6)

1.2 Ansys中的J積分計算策略

沿定義路徑的積分算法早已實現(xiàn)，以2D面積積分為例，其近似為沿Γi的n個線積分之和，乘以環(huán)段wi的寬度，如圖6所示。實際上，wi在這里被設置為一個常數(shù)，這在利用有限元網(wǎng)格計算時是有用的?？梢钥闯?，對裂紋尖端進行了加密細化，離裂紋尖端越遠，寬度隨之增加。Ansys中的積分計算策略認為[17]，最內側和最外側積分環(huán)的寬度必須適合積分區(qū)域內最小和最大的有限元網(wǎng)格。

圖6 Ansys中進行區(qū)域積分的不同積分路徑Fig.6 The different integral paths of integral region in Ansys

(7)

(8)

(9)

(10)

(11)

(12)

在3D情況下，計算不同ξ位置的3個面積積分，乘以加權因子δsj，然后相加。沿ξ3的積分格式類似于一維高斯積分，位置ξj和權重δsj在式(11)(12)中表示。應注意圓柱形積分體積δs的深度，因為它們必須等于或小于沿裂紋前緣的單元尺寸。

1.3 含裂紋RPV有限元模型

本文重點研究RPV內熔融物滯留情況下的斷裂力學行為。在有限元建模中，外部水冷效果由溫度相關傳熱函數(shù)給出，在IVR啟動注水后，RPV外部將被快速淹沒，注水液面上升比較快，因此本文分析中假設RPV外表面被冷水淹沒的情況。在臨界熱流密度下，從核狀沸騰到膜狀沸騰過渡，RPV的外表面最高溫度約為150 ℃[9]。由于單元(Solid 182)不允許直接熱-機耦合，因此必須執(zhí)行順序耦合，計算分解為預定義的時間步長，對于每個時間步，得到一個瞬態(tài)解，并且在隨后的力學計算中采用空間溫度場作為體載荷來計算熱相關場，以及應力場，文獻[18]證明以上有限元計算方案有利于結果收斂。為此，本文利用Ansys計算了3種不同內外徑比(Ro/Ri=1.01，1.10，1.70)RPV上不同半橢圓裂紋排列結構下沿裂紋前緣的3DJ積分。

1.4 半橢圓表面裂紋的3D子模型

建立半橢圓形表面裂紋的3D子模型，位置和形式如圖7所示。由于瞬態(tài)溫度場由熱計算求得，所以只需要得到力學解。這種方法在前面已描述，用于獲得熱機械解，并使用第1.2節(jié)中的3D積分模型和計算策略來計算沿3D裂紋前緣的J積分。圖7示出了帶有裂紋的RPV模擬中的邊界和起始條件，3D軸對稱模型在外側的最大應力位置(也稱為“焊縫”)處包含周向裂紋。對于邊界條件，外表面溫度為150 ℃，而RPV內熔體溫度約為1 327 ℃，與臨界熱流條件相對應。

圖7 RPV含裂紋的三維有限元模型Fig.7 The 3D finite element model of crack-containing RPV

1.5 裂紋的無量綱J積分及驗證

當工作溫度較大時，應考慮材料的延性，采用修正的彈性本構關系作為基準溫度下延性斷裂分析的基礎。在韌性斷裂模型中，損傷起始準則為最大主應力準則，演化規(guī)律以能量釋放率為判據(jù)，彈塑性斷裂力學(EPFM)分析采用指數(shù)軟化模型。在平面應變條件下，能量釋放率可由下式計算：

(13)

根據(jù)文獻[19]給出的K00，所有的歸一化應力強度因子(SIF)均進行歸一化處理，薄、厚球形壓力容器外表面產(chǎn)生的環(huán)向裂紋的SIF與裂紋深度的關系如下:

(14)

(15)

Q為橢圓裂紋的形狀因子，由第二類完整橢圓積分的平方給出，通常近似[20]為：

(16)

圖8 有限元法和API 579標準中解的對比Fig.8 Comparison of J-integral solutions of finite element and API 579

由于幾何構型的對稱性不同，只對含裂紋的球形下封頭進行了分析，赤道面ψ=0°和縱面ψ=90°如圖9中定義。下封頭內表面施加p=2.5 MPa壓力，在裂紋張開面考慮內壓的作用。為了適應裂紋前緣附近的奇異應力場，采用一層20節(jié)點等參塊狀單元折疊成楔狀，在裂紋前緣形成奇異單元，并劃分相應位置，頂部至少有4層網(wǎng)格。此外，模型的其余部分同時包含了20節(jié)點塊狀和10節(jié)點四面體單元。裂紋尖端附近的網(wǎng)格很小，網(wǎng)格密度隨裂紋尖端距離的增大而增大，如圖9所示。利用裂紋面位移外推法提取J積分，并將其計算方案以用戶子程序構建在Ansys中。考慮到RPV上最可能的裂紋位置，假設所有裂紋在赤道面上是相同、等間距和共面的。在有限元計算中，考慮了具有不同a/c的最可能的半橢圓形裂紋，裂紋形式如圖9所示。

圖9 不同橢圓度裂紋前緣區(qū)域的有限元網(wǎng)格Fig.9 The grids of crack front edges of different semi-elliptical cracks

2 結果與討論

2.1 RPV場參數(shù)演化及分布

堆芯熔融嚴重事故條件下，冷卻水最終淹沒RPV，因為水注滿整個堆腔的時間相對較短，本文不考慮水位上升的過程。圖10示出了不同時間沿著壁厚路徑A的溫度分布，可以看出，內壁溫度等于開始時的熔池溫度，并且隨著時間的增加而緩慢下降。實際上，在前幾秒內，靠近內壁的溫度有急劇下降的趨勢，而隨著時間的增加，整個壁厚仍存在較高的溫度梯度。雖然72 h內的內壁溫度略低于10 h時的內壁溫度，但t≥1 h時沿壁厚溫度分布已非常相近。

圖10 不同時間下路徑A沿壁厚的溫度分布Fig.10 Temperature distribution of path A along wall thickness at different time

從RPV外表面(x=149 mm)上看，溫度梯度實際上是隨著時間的增加而上升的，在3 600 s(1 h)左右達到最大值。因此，該溫度梯度導致容器外部出現(xiàn)較大的環(huán)向應力，如圖11所示。

圖11 不同時間路徑A沿壁厚的應力分布情況Fig.11 Stress distribution of path A along wall thickness at different time

圖11揭示了壓縮狀態(tài)應力朝著內壁表面減小，導致其應力承載能力減少，主要是因為熔融物接近或高于RPV材料熔點的高溫作用。此外，最大應力出現(xiàn)在球形封頭和圓柱段的過渡區(qū)，即焊縫位置，這個位置容易出現(xiàn)裂縫，也是最危險的地方。因此，假設半橢圓形裂紋出現(xiàn)在上述位置。

在傳統(tǒng)的IVR概念中，RPV的內表面暴露于溫度T>1 000 ℃的熔池中，而內壓是通過多級泄壓閥釋放完全，因此只要RPV材料沒有被熔穿，即不會發(fā)生所謂的“熱失效”，堆芯的熔融物就被滯留在RPV內。但是實際的事故情況下存在內壓和結構裂紋的影響，這時RPV的破裂風險急劇加大。圖12示出在p=2.5 MPa內壓下含半橢圓形裂紋RPV的應力場和溫度場，由此可知，裂紋尖端處應力集中現(xiàn)象非常明顯，已超過材料的屈服強度，在內壓和熱應力的作用下發(fā)生了整體膨脹和局部外凸，特別是含裂紋區(qū)域存在更大的應力和溫度梯度。

圖12 半橢圓形裂紋周圍的應力及溫度分布Fig.12 Stress and temperature distribution around the semi-elliptical crack

圖13(a)示出不同橢圓度裂紋(a/c=0.3，1.0，1.5)前緣周圍的等效應力分布。結果表明，裂紋尖端的應力集中現(xiàn)象比較明顯，但裂紋尖端的應力集中并不總是發(fā)生在裂紋最深處。實際上，圖13描述的周向裂紋處于1/4對稱RPV模型外側應力最大處，應力最大值約680 MPa。從圖13還可以看出，3種橢圓度裂紋的整體應力水平不一致，以裂紋最深處應力水平看，呈現(xiàn)a/c=0.3的應力水平大于a/c=1.0的，同時a/c=1.0的應力水平大于a/c=1.5的。從裂紋剖切面看，應力分布明顯不對稱，如圖13(b)所示。值得注意的是，不同橢圓度裂紋最大應力處也不同，比如針對a/c=0.3橢圓度裂紋，最大應力處出現(xiàn)在靠近中間深度，而對于a/c=1.0，1.5的裂紋，最大應力處出現(xiàn)在靠近RPV外表面的裂紋尖端。不同橢圓度裂紋對應不同的裂紋前緣應力分布，a/c=1.0橢圓度對應的裂紋前緣應力較小，應力集中現(xiàn)象相對不明顯，而a/c=0.3橢圓度的RPV出現(xiàn)最嚴重的應力集中現(xiàn)象。

(a)水平視角

2.2 裂紋參數(shù)對無量綱J積分的影響

為了確定不同參數(shù)對裂紋前緣J積分分布的影響，考慮了半橢圓裂紋的3種幾何形式，橢圓度分別為a/c=0.3,1.0,1.5,不同裂紋深度a/t=0.05～0.55，以及不同裂紋密率陣列δ=0～0.9情況下的裂紋前緣J積分分布，分析評估了JI/Jo分布情況，其中Jo為a=15 mm情況下的J積分結果。

典型球形下封頭Ro/Ri=1.1，當裂紋密率δ=0～0.9時，在相對淺裂紋a/t=0.3(見圖14(a))和深裂紋a/t=0.6(見圖14(b))前緣的J/Jo隨參數(shù)角ψ的變化如圖14所示。由圖14(a)(b)的比較可以看出，沿裂紋前緣的JI/Jo分布具有不同的模式，a/t=0.3淺裂紋的J/Jo分布的模式與a/t=0.6深裂紋的幾乎相反。對于淺裂紋情況，JI/Jo的最大值出現(xiàn)在ψ=90°的最深裂紋處，如圖14(a)所示；而對于深裂紋情況，JI/Jo的最大值發(fā)生在尖端ψ=0°附近，如圖14(b)所示。對于淺裂紋或深裂紋，沿整個裂紋前緣的JI/Jo值隨裂紋密率δ的增加而增大。在共面周向裂紋中，沿整個裂紋前緣的JI/Jo隨裂紋密率δ的增加而增大，反映了陣列裂紋的增加導致RPV承載力的下降，裂紋往前擴展加劇。相鄰裂紋之間的主要相互作用發(fā)生在ψ=0°的尖端，從而導致JI/Jo梯度和Jmax增大。

(a)a/t=0.3

為了評估下封頭裂紋密率對裂紋前緣最大J積分的影響，對a/c=0.3的橢圓度和Ro/Ri=1.1的幾何形狀，δ=0.9，0.7，0的3個共面裂紋陣列的Jmax/Jo與裂紋深度的關系分析如圖15所示。盡管裂紋密率δ有所不同，但Jmax/Jo與a/t的之間關系趨勢相似，如圖15所示，隨著δ的增大，Jmax/Jo隨之增大。此外，隨著裂紋深度a/t的增加，Jmax/Jo單調增加，Jmax/Jo與裂紋深度a/t呈分段線性關系。由圖15可以看出，Jmax/Jo與a/t之間關系曲線可分為兩個階段，即初期快速增加和后期穩(wěn)定增加。

圖15 橢圓度a/c=0.3下不同裂紋密率δ對Jmax/Jo變化影響Fig.15 Effect of crack density δ on Jmax/Jo for crack a/c=0.3

與橢圓度a/c=0.3裂紋J分布相比,a/c=1.0的裂紋最大JI的位置在尖端ψ=0°與裂紋最深點ψ=90°之間，并靠近ψ=0°的位置，JI隨裂紋深度的增加而呈現(xiàn)不同分布，如圖16所示。對于較淺裂紋和較深裂紋，JI值也隨裂紋密率的增加而加快上升，沿整個裂紋前緣的JI/Jo分布隨裂紋密率δ的增加而增大。在a/c=1.0的裂紋橢圓度下，a/t=0.3和a/t=0.6深度的JI/Jo沿裂紋前緣分布相似。由圖16可以看出，當裂紋深度變大時，最大JI所在位置更接近最尖端點ψ=0°，而且最大JI值更高。與最大JI位置不同，對于淺裂紋和深裂紋，最小JI都出現(xiàn)于最內側點ψ=90°(也稱為“最深裂紋深度”)。JI/Jo隨著裂紋深度a/t的增加而增加，但隨著裂紋深度的增加，裂紋密率δ引起的JI/Jo數(shù)據(jù)的離散更為明顯，如圖16(b)所示。由圖16(b)還可以看出，當a/t=0.6和δ=0.9時，最大JI是Jo的7.0倍，而a/t=0.6和δ=0時，最大JI是Jo的3.2倍。從圖17可以看出，當裂紋形狀接近圓形(a/c=1.0)時，裂紋密率δ對Jmax/Jo的影響非常敏感，而且影響也很大，隨著裂紋深度a/t增加，Jmax/Jo單調增加。

(a)a/t=0.3

圖17 橢圓度a/c=1.0下不同裂紋密率δ對Jmax/Jo變化影響Fig.17 Effect of crack density δ on Jmax/Jo for crack a/c=1.0

相對深度分別為a/t=0.3和0.6，對于橢圓度a/c=1.5裂紋，最大JI/Jo總是位于裂紋尖端ψ=0°附近，如圖18所示。對比圖18(a)和(b)可以看出，對于相對較淺裂紋a/t=0.3，最大JI/Jo附近的集中分布現(xiàn)象更為顯著，如圖18(a)所示。實際上，最大JI/Jo一致出現(xiàn)在裂紋尖端ψ=0°附近，并朝著尖端略微減小，如圖18(b)所示。

圖18和圖16相比，a/c=1.5對應的沿裂紋前緣的JI/Jo分布類似于a/c=1.0(圓形裂紋輪廓)的JI/Jo分布。盡管如此，圖18揭示了由于較高的裂紋密率δ導致整體上JI/Jo相對較大。此外，隨著裂紋密率δ的增加，尖端點附近JI/Jo的下降變得更加顯著。比較圖18(a)和(b)可以看出，隨著裂紋深度a/t的增加，JI/Jo的值雖然變大，但是JI/Jo峰值變得更加平坦。由圖18可以看出，無量綱J積分的數(shù)據(jù)離散度隨橢圓表面裂紋深度a/t的增加而增大，此外，還發(fā)現(xiàn)JI/Jo數(shù)據(jù)離散度隨裂紋密率δ的增加而增大。J積分的數(shù)據(jù)離散主要是因為RPV壁厚上的溫度分布不均勻，加上沿壁厚方向一次應力、二次應力和峰值應力的共同作用。

(a)a/t=0.3

橢圓度a/c=1.5下不同裂紋密率δ對Jmax/Jo隨a/t變化的影響如圖19所示，可以看出，裂紋密率δ對Jmax/Jo-a/t的關系有較大的影響?？傮w上看，Jmax/Jo-a/t關系曲線還是呈現(xiàn)分段線性和單調增加的特征，另外，隨著裂紋密率δ的增加，Jmax/Jo表現(xiàn)為先緩慢增加(0<δ<0.7)、后快速增加(0.7<δ<0.9)。與圖15，17不同的是，δ=0.7和δ=0.9對應Jmax/Jo-a/t曲線之間的差值更小，即a/c=1.5橢圓度裂紋的Jmax/Jo絕對值相對a/c=0.3，1.0橢圓度裂紋的Jmax/Jo較小，意味著a/c=1.5橢圓度裂紋對裂紋密率δ的敏感性較低，同樣載荷條件下裂紋更不容易擴展。綜合來看，同樣的裂紋深度a/t下，a/c=1.5橢圓度裂紋在不同的裂紋密率δ下Jmax/Jo的分散性最小，即小于a/c=0.3，1.0橢圓度裂紋Jmax/Jo的分散性。

圖19 橢圓度a/c=1.5下不同裂紋密率δ對Jmax/Jo隨a/t變化的影響Fig.19 Effect of crack density δ on Jmax/Jo with changing of a/t crack a/c=1.5

3 結論

(1)利用熱-機耦合下斷裂力學方法，研究了堆芯熔毀事故下反應堆壓力容器(RPV)下封頭的斷裂行為。計算了含半橢圓形裂紋的3D裂紋的J積分，考慮了裂紋橢圓度分別為a/c=0.3，1.0，1.5情況。結果顯示，對于不同的半橢圓形裂紋，裂紋前緣的J積分分布顯著不同，尤其是a/c=0.3的裂紋，其J積分分布與其他兩種相反。對于a/c=0.3的半橢圓形裂紋，隨著裂紋深度的增加，J積分分布會發(fā)生本質性的變化。對于較淺裂紋(a/t=0.3)，J積分最大值出現(xiàn)在最大裂紋深度處，而對于較深裂紋(a/t=0.6)，J積分最大值出現(xiàn)在ψ=0°的裂紋尖端附近。

(2)裂紋密率δ和裂紋橢圓度a/c是影響RPV壓力邊界完整性的兩個重要因素。隨著δ的減小，J積分峰值明顯減弱。RPV周圍材料對裂紋的穩(wěn)定作用更為顯著。一般，δ=0的J積分約為δ=0.9的J積分的2～3倍。

(3)如果不考慮RPV壁厚上的高溫度梯度，勢必會錯誤估算J積分分布和峰值大小，可能會較大地高估和低估J積分值，說明在嚴重事故條件下對RPV進行斷裂力學評價時，熱-機邊界條件精度對RPV安全設計的重要性。