大破口失水事故分析方法GSM的開發與應用

曹志偉，林支康，王 婷，梁 任，鮑 杰，盧向暉

(1.中廣核研究院有限公司，廣東 深圳 518031；2.生態環境部 核與輻射安全中心，北京 100082)

大破口失水事故(LB LOCA, large break loss of coolant accident)是核電廠最重要的設計基準事故之一，是驗證應急堆芯冷卻系統(ECCS)能力并直接限制核電廠發電能力的極限事故，LB LOCA裕量也是最重要的核電廠安全裕量指標之一。由于LB LOCA現象復雜，后果嚴重，在核電廠安全分析中處于非常重要的地位。LB LOCA分析程序、方法及分析均是核電廠建造和裝料許可證等申請的關鍵內容。目前擁有的LB LOCA分析方法是20世紀90年代開發的較為保守的分析方法，保守的分析方法掩蓋了核電廠的真實安全裕量，如18個月換料模式下的CPR1000核電廠，最低LOCA裕量只有1%左右。保守分析方法無法真實反映核電廠設計的安全指標和先進性，也極大地限制了核電廠經濟性的提升。因此，為提升核電廠安全裕量，本文基于現有的最佳估算熱工水力分析程序CATHARE GB及相關保守方法，開發基于統計理論的先進LOCA分析方法，并初步應用于CPR1000核電廠。

1 確定論現實方法介紹

20世紀90年代，法國電力公司和法瑪通公司基于其多年的核電廠設計和運營經驗，將LB LOCA分析方法逐步從保守的、滿足美國聯邦法規10CFR50附錄K[1]的方法論轉向了基于最佳估算程序的分析方法——確定論現實方法(DRM, deterministic realistic method)[2]。程序計算結果的不確定性來源于兩方面：初始邊界條件和軟件物理模型的不確定性。分析所用的計算程序是現實的，即程序必須對所有瞬態物理現象進行最佳估算。DRM開發基于CATHARE 2 V1.3L程序，針對LB LOCA現象提供了專用的模型。這些模型均與現有實驗進行對比，得到程序計算模型的不確定性。基于現實的電廠模型，通過對程序計算模型的不確定性、電廠技術參數及初始狀態的不確定性進行抽樣計算，得到每個目標參數(如燃料包殼峰值溫度，PCT)在95%置信水平下、95%概率下不會被目標參數值超過的保守值(即雙95%值)。隨后定義出程序的保守模型(DRM懲罰模型)，該模型計算結果可包絡目標參數的雙95%值，因此通過DRM得到的目標參數結果高于統計方法得到的95%置信水平的結果。所構造的懲罰模型不僅要求保守性最小，還應保留程序對真實瞬態工況的模擬能力。DRM的開發主要包括4個方面：1) 評估模型的現實性；2) 量化不確定性；3) 引入懲罰模型；4) 保守性評價。

在DRM中為了滿足計算結果的保守性，引入了程序懲罰模型。為了盡可能地既不影響瞬態過程，又能同時包絡PCT95/95值，最終選擇對幾個參數進行懲罰，稱為DRM懲罰模型[1]。

2 確定論現實方法開發

在LOCA保守評價模型中，影響LOCA評估的相關模型和電廠狀態依據美國聯邦法規10CFR50附錄K[1]的要求均須采用保守方式處理。在LOCA最佳估算評價模型[3-4]中，程序物理模型和電廠狀態參數的不確定性采用統計方法處理，以綜合呈現出主要評估參數的不確定性范圍。影響LOCA最佳估算不確定性的因素分為兩大類：一類是程序模型的不確定性，另一類是電廠狀態參數的不確定性。

中廣核確定論統計方法(GSM, CGN deterministic statistic methodology)為一種介于保守評價模型和最佳估算評價模型兩類方法之間的折衷式LOCA分析方法。在該方法中，重要程序模型部分以經實驗驗證的保守方法處理；而電廠重要狀態參數的處理，則遵循CSAU的執行程序，采用統計方法量化電廠重要狀態參數不確定性的整體效應。對于無法實行不確定性分析的電廠重要參數則采用包絡值，以保守處理該參數的不確定性影響。

GSM是基于現有的最佳估算熱工水力分析程序CATHARE GB以及DRM分析相關懲罰模型假設而開發的LOCA分析方法，在GSM中，考慮的不確定性主要有3類：1) CATHARE GB程序模型相關不確定性；2) 電廠模型相關假設；3) 電廠狀態相關參數的不確定性。

2.1 程序模型相關不確定性

對于CATHARE GB程序模型相關不確定性，采用了DRM懲罰模型，相關參數和取值列于表1。表1中：HTCBerenson為噴放階段堆芯傳熱系數；HTCReflooding為再淹沒階段堆芯傳熱系數；Tmsfb為最小膜態沸騰轉變溫度；Qle為氣液兩相間的冷凝系數。

表1 懲罰模型相關參數和取值Table 1 Relevant parameter and value of penalization model

2.2 電廠模型相關假設

電廠模型相關的不確定性包括：考慮單一故障準則；破口發生時刻考慮喪失廠外電源，主泵停運；事故發生時的燃料循環，采用保守的燃耗數據，通過敏感性分析得到；一回路每個環路的流量為熱工水力設計流量；考慮最大的堆芯旁流量，使通過堆芯的冷卻劑量最少；一回路質量計算考慮蒸汽發生器堵管率；低安全殼壓力使得再淹沒進程減緩并降低再淹沒期間堆芯內熱傳遞系數，對PCT更加不利，因此在計算中采用正常運行工況下最小的初始安全殼溫度和壓力，并選擇安全殼最大濕度。破口尺寸和類型采用統計抽樣處理，其中，破口類型為等概率分布，破口面積的抽樣范圍為[0.241 6A，2.0A]，A為冷管段的橫截面積。

2.3 電廠狀態參數的不確定性

電廠狀態參數的不確定性處理重點是要篩選出重要電廠狀態參數并對其進行不確定性量化處理。

GSM對于電廠狀態相關參數采用了統計抽樣處理，針對CPR1000核電廠，采用抽樣統計處理的電廠狀態相關參數列于表2。

表2 重要電廠參數及其不確定性范圍和分布形式Table 2 Plant status parameter and it’s uncertainty range and distribution

2.4 重要參數的隨機抽樣處理

用統計方法量化計算誤差時，必須針對相關不確定性因素進行隨機取樣，產生足夠數目的隨機樣本。隨機樣本數目由所采用的統計方法決定。一般而言，統計可分為非參數方法和參數方法兩大類。非參數方法的統計分析不需要知道分析結果誤差概率密度分布類型，而參數方法的統計分析必須假設分析結果誤差概率密度分布特性，如函數種類、均值及標準方差等[5]。GSM采用簡單隨機抽樣，對電廠參數進行統一抽樣，得到59組統計抽樣組合工況并計算得到59個PCT值，供上述兩種統計方法分析使用，計算PCT95/95統計上限值。

2.5 重要目標參數不確定性量化統計分析

采用考慮DRM保守懲罰的CATHARE GB程序模型對59組工況計算得到59組PCT，并對此進行統計處理以量化在考慮電廠狀態參數不確定性后的整體計算結果的誤差范圍，并進一步計算或估算出PCT95/95統計上限值。采用的統計處理方法有兩種，即正態分布單側置信上限方法(參數統計法)和順序統計方法(非參數統計法)。

1) 正態分布單側置信上限方法[6]

在此統計分析方法中先假設59組工況的分析結果滿足正態分布，并經適當方法檢驗確認，如χ2擬合檢驗法。如果滿足特定函數分布，則可由樣本均值μs及樣本標準方差σs在一定的置信度95%下估算總體均值μp及標準方差σp，并進一步估算上限值，如PCT95/95之值。樣本均值和樣本標準方差分別為：

其中：xi為樣本變量參數值；n為樣本數量。

參數統計分析可分為3個步驟，即正態分布檢驗、由樣本均值與標準方差估算總體均值與標準方差及由分析結果計算95/95統計上限值。

(1) 正態分布檢驗

所采用的正態分布檢驗方法是χ2擬合檢驗法，該方法有兩項適用要求：只適合大樣本情形，一般要求樣本容量N≥50；每個區間的樣本數(NPi)不應太小，其中Pi為子區間i對應的概率，一般要求不小于5，否則應適當合并鄰近區間，使得NPi≥5，同時區間數k一般應在4～20之間。當假設分布為正態分布N(μ，σ)時，真實的總體μ、σ仍是未知。根據極大似然估計法可由樣本μs和σs合理估算。

將樣本數NPi依大小分為數個區間，定義正態分布的分歧度υ，量化各區間預期樣本數和實際樣本數的差異：

(2) 由樣本均值與樣本標準方差估算總體均值與標準方差

以PCT上限包絡值的計算為例，總體均值μp在1-α置信度下單側置信上限為：

Pμp≤μs+tα(n-1)σs/n{}=1-α

其中，tα為學生氏分布統計參數。

同樣地，總體標準方差σp在1-α置信度下單側置信上限為：

(3) 計算95/95統計上限值

上述相關計算式中在一定置信水平或置信度1-α下，μp和σp上限包絡值可由樣本均值和標準方差直接計算。當μp和σp上限包絡值在95%置信度下計算出后，目標參數的95/95統計上限值Y95/95可直接由下式計算：

Y95/95=μp,95%+1.645σp,95%

其中，μp,95%和σp,95%分別為在95%置信度下μp和σp的上限包絡值。

2) 非參數統計法

采用Wilk’s非參數統計法[7]不需先行確定樣本所服從的概率分布，而是在有限的隨機抽樣樣本中，如n=59、93或124組，保守估算可能的95/95統計上限值：Y95/95≈Y1st(59)、Y95/95≈Y2nd(93)或Y95/95≈Y3rd(124)。

3) 最終目標參數值的確定

在確定最終PCT值時，如果正態分布單側置信上限方法和非參數統計法均能計算或估算出PCT95/95值，用較保守者作為最終PCT值，即：

PCT=max(PCT95/95,PCT非參數統計)

3 GSM應用

將GSM應用于CPR1000核電廠某燃料管理方案研究項目中，其分析結果與傳統保守的DRM進行比較，以合理挖掘LB LOCA的裕量。不確定性分析的重要電廠參數范圍和分布列于表2。經統一隨機抽樣59組樣本組合工況，并使用CATHARE GB程序模型計算得到59組樣本對應的PCT值。GSM應用的流程如圖1所示。

3.1 非參數統計處理

根據Wilk’s公式γ=1-βn，當概率水平(γ)、置信水平(β)均設定為95%時，n即為59。此時59組抽樣結果的最大值即為95/95統計上限值的保守估算結果。采用Wilk’s非參數統計方法計算得到PCT95/95統計上限值為1 136.63 ℃。

3.2 參數統計處理

參數統計分析分為3個步驟：1) 正態分布檢驗；2) 由樣本均值與標準方差估算總體均值與標準方差；3) 由分析結果計算95/95統計上限值。

1) 正態分布檢驗

圖1 GSM流程Fig.1 Procedure of GSM

首先樣本均值和標準方差的計算如下：

2) 由μs和σs估算總體均值與標準方差

μp的95%單側置信區間為(-∞，899.244]，σp的95%單側置信區間為(-∞，94.514]。

由于滿足正態分布，所以95/95的PCT值為：

PCT95/95=μp,95%+1.645σp,95%=

899.244+1.645×94.514=1 054.72 ℃

表3 擬合優度檢驗分歧度表Table 3 Data for goodness-of-fit test

3) 計算95/95統計上限值

由Wilk’s非參數統計法得到PCT1st=1 136.63 ℃，最終計算得到PCT95/95的值為：

PCTGSM=max(PCT1st,PCT95/95)=1 136.63 ℃

GSM應用于CPR1000核電廠LB LOCA分析結果與傳統DRM分析結果的比較如圖2所示。由圖2可見，采用GSM計算得到的LB LOCA分析結果比傳統DRM分析結果低117 ℃，可提高約9%的裕量。

圖2 GSM與DRM分析結果對比Fig.2 Comparison of analysis result of GSM and DRM

4 小結

GSM是介于保守評價模型和最佳估算評價模型兩類方法之間的折衷式LOCA分析方法。在該方法中，程序物理模型仍以DRM懲罰模型進行保守方法處理，對核電廠重要狀態參數采用統計方法量化各重要電廠狀態參數不確定性范圍和分布，并對統計抽樣計算得到的目標參數(PCT)進行統計處理以得到PCT95/95。GSM應用于CPR1000核電廠LB LOCA分析可挖掘約9%的LOCA裕量。