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

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

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

1 確定論現(xiàn)實(shí)方法介紹

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

在DRM中為了滿足計(jì)算結(jié)果的保守性，引入了程序懲罰模型。為了盡可能地既不影響瞬態(tài)過程，又能同時(shí)包絡(luò)PCT95/95值，最終選擇對幾個(gè)參數(shù)進(jìn)行懲罰，稱為DRM懲罰模型[1]。

2 確定論現(xiàn)實(shí)方法開發(fā)

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

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

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

2.1 程序模型相關(guān)不確定性

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

表1 懲罰模型相關(guān)參數(shù)和取值Table 1 Relevant parameter and value of penalization model

2.2 電廠模型相關(guān)假設(shè)

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

2.3 電廠狀態(tài)參數(shù)的不確定性

電廠狀態(tài)參數(shù)的不確定性處理重點(diǎn)是要篩選出重要電廠狀態(tài)參數(shù)并對其進(jìn)行不確定性量化處理。

GSM對于電廠狀態(tài)相關(guān)參數(shù)采用了統(tǒng)計(jì)抽樣處理，針對CPR1000核電廠，采用抽樣統(tǒng)計(jì)處理的電廠狀態(tài)相關(guān)參數(shù)列于表2。

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

2.4 重要參數(shù)的隨機(jī)抽樣處理

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

2.5 重要目標(biāo)參數(shù)不確定性量化統(tǒng)計(jì)分析

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

1) 正態(tài)分布單側(cè)置信上限方法[6]

在此統(tǒng)計(jì)分析方法中先假設(shè)59組工況的分析結(jié)果滿足正態(tài)分布，并經(jīng)適當(dāng)方法檢驗(yàn)確認(rèn)，如χ2擬合檢驗(yàn)法。如果滿足特定函數(shù)分布，則可由樣本均值μs及樣本標(biāo)準(zhǔn)方差σs在一定的置信度95%下估算總體均值μp及標(biāo)準(zhǔn)方差σp，并進(jìn)一步估算上限值，如PCT95/95之值。樣本均值和樣本標(biāo)準(zhǔn)方差分別為：

其中：xi為樣本變量參數(shù)值；n為樣本數(shù)量。

參數(shù)統(tǒng)計(jì)分析可分為3個(gè)步驟，即正態(tài)分布檢驗(yàn)、由樣本均值與標(biāo)準(zhǔn)方差估算總體均值與標(biāo)準(zhǔn)方差及由分析結(jié)果計(jì)算95/95統(tǒng)計(jì)上限值。

(1) 正態(tài)分布檢驗(yàn)

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

將樣本數(shù)NPi依大小分為數(shù)個(gè)區(qū)間，定義正態(tài)分布的分歧度υ，量化各區(qū)間預(yù)期樣本數(shù)和實(shí)際樣本數(shù)的差異：

(2) 由樣本均值與樣本標(biāo)準(zhǔn)方差估算總體均值與標(biāo)準(zhǔn)方差

以PCT上限包絡(luò)值的計(jì)算為例，總體均值μp在1-α置信度下單側(cè)置信上限為：

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

其中，tα為學(xué)生氏分布統(tǒng)計(jì)參數(shù)。

同樣地，總體標(biāo)準(zhǔn)方差σp在1-α置信度下單側(cè)置信上限為：

(3) 計(jì)算95/95統(tǒng)計(jì)上限值

上述相關(guān)計(jì)算式中在一定置信水平或置信度1-α下，μp和σp上限包絡(luò)值可由樣本均值和標(biāo)準(zhǔn)方差直接計(jì)算。當(dāng)μp和σp上限包絡(luò)值在95%置信度下計(jì)算出后，目標(biāo)參數(shù)的95/95統(tǒng)計(jì)上限值Y95/95可直接由下式計(jì)算：

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

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

2) 非參數(shù)統(tǒng)計(jì)法

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

3) 最終目標(biāo)參數(shù)值的確定

在確定最終PCT值時(shí)，如果正態(tài)分布單側(cè)置信上限方法和非參數(shù)統(tǒng)計(jì)法均能計(jì)算或估算出PCT95/95值，用較保守者作為最終PCT值，即：

PCT=max(PCT95/95,PCT非參數(shù)統(tǒng)計(jì))

3 GSM應(yīng)用

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

3.1 非參數(shù)統(tǒng)計(jì)處理

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

3.2 參數(shù)統(tǒng)計(jì)處理

參數(shù)統(tǒng)計(jì)分析分為3個(gè)步驟：1) 正態(tài)分布檢驗(yàn)；2) 由樣本均值與標(biāo)準(zhǔn)方差估算總體均值與標(biāo)準(zhǔn)方差；3) 由分析結(jié)果計(jì)算95/95統(tǒng)計(jì)上限值。

1) 正態(tài)分布檢驗(yàn)

圖1 GSM流程Fig.1 Procedure of GSM

首先樣本均值和標(biāo)準(zhǔn)方差的計(jì)算如下：

2) 由μs和σs估算總體均值與標(biāo)準(zhǔn)方差

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

由于滿足正態(tài)分布，所以95/95的PCT值為：

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

899.244+1.645×94.514=1 054.72 ℃

表3 擬合優(yōu)度檢驗(yàn)分歧度表Table 3 Data for goodness-of-fit test

3) 計(jì)算95/95統(tǒng)計(jì)上限值

由Wilk’s非參數(shù)統(tǒng)計(jì)法得到PCT1st=1 136.63 ℃，最終計(jì)算得到PCT95/95的值為：

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

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

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

4 小結(jié)

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