船用壓水堆中破口失水事故不確定性及敏感性分析*

王海峰 陳玉清

（海軍工程大學(xué) 武漢 430033）

1 引言

早期失水事故（LOCA）分析采用保守的程序和各種偏保守的假設(shè)來確保安全，但過于保守的假設(shè)不一定能得出保守的結(jié)果，且不利于建立合理的事故處置規(guī)程。1988 年，NRC 對(duì)10CFR50.46 及其附錄K 進(jìn)行了修訂，允許采用最佳估算加不確定性（BEPU）分析方法進(jìn)行LOCA 分析，但需確保驗(yàn)收準(zhǔn)則規(guī)定的限值有很高的概率不被超過。1989年，NRC發(fā)布了管理導(dǎo)則RG1.157［1］，對(duì)最佳估算程序的計(jì)算模型、經(jīng)驗(yàn)關(guān)系、取用數(shù)據(jù)、評(píng)價(jià)方法做出明確的規(guī)定，并開發(fā)出程序比例、適用性和不確定性評(píng)價(jià)方法CSAU（Code Scaling，Applicability，and Uncertainty）［2］。

針對(duì)壓水堆失水事故，曾未等對(duì)壓水堆失水事故實(shí)驗(yàn)（LOFT）的FP-LP-2 工況進(jìn)行了模擬計(jì)算，并采用德國(guó)反應(yīng)堆安全研究所（GRS）不確定性分析方法對(duì)計(jì)算結(jié)果進(jìn)行不確定性和敏感性分析，給出了關(guān)鍵輸入?yún)?shù)95%置信度的不確定性包絡(luò)帶，確認(rèn)了GRS 方法的有效性［3］。陳玉清等對(duì)船用壓水堆小破口失水事故進(jìn)行了不確定性和敏感性分析，驗(yàn)證了事故進(jìn)程中反應(yīng)堆可以實(shí)現(xiàn)堆芯的自我冷卻，確定了對(duì)包殼峰值溫度（PCT）較為敏感的輸入?yún)?shù)［4］。楊曉敏等采用最佳估算加不確定性分析的方法對(duì)小型壓水堆小破口失水事故（SB-LOCA）開展了瞬態(tài)特性分析，確定了對(duì)PCT 影響較大的關(guān)鍵不確定性源系數(shù)［5］。

相比于SB-LOCA，MB-LOCA 進(jìn)程中冷卻劑裝量減少更快、事故處置時(shí)間更短、對(duì)反應(yīng)堆安全威脅更大。目前針對(duì)船用壓水堆MB-LOCA 分析研究較少，本文基于輸入不確定性傳播的GRS 方法，采用RELAP5/MOD3.2最佳估算程序建立全系統(tǒng)熱工水力計(jì)算模型，模擬48mm 中破口失水事故瞬態(tài)進(jìn)程，對(duì)運(yùn)算結(jié)果開展不確定性量化計(jì)算和敏感性分析，提出事故緩解優(yōu)化策略，對(duì)優(yōu)化事故處置具有積極作用。

2 不確定性分析方法

目前國(guó)際上普遍采用BEPU 分析方法進(jìn)行LOCA 研究，方法主要分為基于輸入的不確定性傳播評(píng)估方法（如CSAU，GRS，ASTRUM）與基于輸出的不確定性傳播評(píng)估方法（如UMAE）兩類［6］，其中，基于輸入的不確定性傳播評(píng)估方法的研究和應(yīng)用相對(duì)更加廣泛。GRS 評(píng)價(jià)方法由德國(guó)反應(yīng)堆安全研究所開發(fā)，以最佳估算程序ATHLET 作為工具，考慮程序模型參數(shù)、初始和邊界條件參數(shù)等輸入?yún)?shù)的不確定性，采用基于Wilks 公式的非參數(shù)統(tǒng)計(jì)方法，進(jìn)行不確定性的傳遞計(jì)算［7］，與CSAU 方法相比，GRS方法的優(yōu)勢(shì)有：1）輸入?yún)?shù)均用取值范圍與概率分布進(jìn)行描述；2）不確定性分析所要計(jì)算的次數(shù)與輸入?yún)?shù)個(gè)數(shù)無關(guān)，只取決于Wilks公式，且次數(shù)遠(yuǎn)小于CSAU 方法采用的響應(yīng)面法；3）通過對(duì)計(jì)算結(jié)果進(jìn)行敏感性分析，得到各輸入?yún)?shù)對(duì)目標(biāo)參數(shù)的影響大小［8］。圖1所示為GRS評(píng)價(jià)方法流程。

圖1 GRS評(píng)價(jià)方法流程

2.1 關(guān)鍵輸入?yún)?shù)及其分布

在壓水堆LOCA 分析中，一般選取PCT 作為評(píng)價(jià)標(biāo)準(zhǔn)。由于最佳估算的LOCA 分析中涉及的參數(shù)較多，且存在大量的不確定性，難以同時(shí)對(duì)所有參數(shù)進(jìn)行不確定性分析，需要基于一定的方法篩選出對(duì)目標(biāo)參數(shù)影響較大的重要輸入?yún)?shù)。目前核電廠壓水堆LOCA 不確定性分析研究較為成熟，主要依據(jù)現(xiàn)象識(shí)別和排序表（PIRT）對(duì)事故現(xiàn)象和參數(shù)進(jìn)行分析［9］，而針對(duì)船用壓水堆MB-LOCA 的不確定性分析應(yīng)用較少，缺少相關(guān)試驗(yàn)數(shù)據(jù)及詳細(xì)信息，因此不確定性輸入?yún)?shù)的范圍及概率分布需要參考電廠堆的研究成果，結(jié)合船用壓水堆安全分析報(bào)告、相關(guān)文獻(xiàn)信息和工程經(jīng)驗(yàn)判斷來確定。

根據(jù)上述對(duì)MB-LOCA 事故進(jìn)程的分析，為保證事故發(fā)生后堆芯的淹沒，綜合考慮船用壓水堆安全注射系統(tǒng)特性和初始及邊界條件對(duì)目標(biāo)參數(shù)的影響，選取7 個(gè)關(guān)鍵不確定性參數(shù)開展不確定性分析，確定其范圍和概率密度分布后對(duì)部分參數(shù)進(jìn)行歸一化處理，結(jié)果如表1 所示。對(duì)于不易量化的不確定性輸入?yún)?shù)，假設(shè)其在相應(yīng)范圍內(nèi)服從均勻分布，且通常認(rèn)為各輸入?yún)?shù)之間相互獨(dú)立。

表1 不確定性輸入?yún)?shù)及概率分布

2.2 輸入?yún)?shù)抽樣

不確定性分析中常采用的抽樣方法有簡(jiǎn)單隨機(jī)抽樣（SRS）和拉丁超立方抽樣（LHS）兩種，LHS方法是SRS方法的改良，其將輸入?yún)?shù)范圍分為等概率且互不重疊的n 個(gè)區(qū)間，并在各區(qū)間中運(yùn)用簡(jiǎn)單抽樣抽取樣本，隨機(jī)組合后再進(jìn)行不確定性傳遞計(jì)算。LHS 比SRS 樣本分布更加均勻，更能體現(xiàn)隨機(jī)變量的分布特征，因此LHS在復(fù)雜系統(tǒng)的不確定性計(jì)算中效率更高。目前求解樣本容量的方法主要有參數(shù)抽樣法和非參數(shù)抽樣法，前者一般要求總體的分布函數(shù)已知，而后者對(duì)總體分布沒有明確要求。LOCA 不確定性分析中，無法統(tǒng)一限定輸入?yún)?shù)的數(shù)量，且程序的輸出結(jié)果（如PCT）一般為未知分布的隨機(jī)變量，故本文采用非參數(shù)統(tǒng)計(jì)中的Wilks方法［10］。根據(jù)Wilks公式可得，單個(gè)輸出變量的單側(cè)容許限值計(jì)算公式為

雙側(cè)容許限值計(jì)算公式為

式中，β為最大計(jì)算結(jié)果不超過限值的置信水平，γ為容許限值的概率份額，N 為計(jì)算次數(shù)。為滿足安全分析“95/95準(zhǔn)則”，當(dāng)β=γ=0.95時(shí)，對(duì)于單側(cè)容忍間N=59，即至少需要成功進(jìn)行59 次抽樣計(jì)算，而對(duì)于雙側(cè)容忍區(qū)間N=93，則至少需要成功進(jìn)行93次計(jì)算。

2.3 模型調(diào)試與計(jì)算

程序建模的范圍主要包括反應(yīng)堆堆芯、蒸汽發(fā)生器、穩(wěn)壓器、泵、閥、管道以及相關(guān)輔助系統(tǒng)、專設(shè)安全系統(tǒng)及控制系統(tǒng)等。進(jìn)行瞬態(tài)計(jì)算前，需要對(duì)所建立的模型進(jìn)行穩(wěn)態(tài)調(diào)試，確保主要輸入?yún)?shù)的計(jì)算值同設(shè)計(jì)值一致，然后采用RELAP5/MOD3.2對(duì)所建模型額定工況條件下MBLOCA 進(jìn)行瞬態(tài)模擬，驗(yàn)證建立的破口模型能夠模擬選定的事故響應(yīng)進(jìn)程。

為滿足不確定性分析計(jì)算次數(shù)的要求，考慮到程序在計(jì)算過程中可能會(huì)出現(xiàn)錯(cuò)誤而終止，本文采用LSH 方法抽取65 組關(guān)鍵輸入?yún)?shù)組合樣本，依次修改輸入卡參數(shù)后運(yùn)行RELAP5 程序模擬事故進(jìn)程，最后對(duì)運(yùn)行結(jié)果開展不確定性和敏感性分析。為提高分析效率，利用Matlab軟件開發(fā)批量處理程序，建立一個(gè)完整有效的失水事故不確定性分析程序，具體流程如圖2所示。

圖2 不確定性分析程序流程圖

3 分析與討論

3.1 不確定性分析

NRC 發(fā)布的輕水堆LOCA 分析驗(yàn)收準(zhǔn)則要求PCT 的最高計(jì)算值不能超過1477K［11］，可將其作為溫度歸一化準(zhǔn)則。程序批量運(yùn)算后得到65 組燃料包殼表面歸一化溫度隨時(shí)間的變化曲線如圖3 所示，反應(yīng)堆堆芯歸一化水位隨時(shí)間的變化曲線如圖4 所示，48mm 中破口失水事故中不同組合抽樣運(yùn)算結(jié)果的燃料包殼表面歸一化峰值溫度分布如圖5所示。

圖3 65組工況下燃料包殼表面溫度

圖4 65組工況下堆芯水位

圖5 48mm中破口失水事故中PCT散點(diǎn)圖

圖3 中各組工況下燃料包殼表面歸一化溫度隨時(shí)間的變化規(guī)律基本一致，但不同的輸入?yún)?shù)組合導(dǎo)致PCT 的高低呈現(xiàn)不確定性差異，圖4 中堆芯歸一化水位隨時(shí)間的變化趨勢(shì)與圖3 相對(duì)應(yīng)，堆芯水位越低，燃料包殼表面溫度越高，且堆芯水位的最低點(diǎn)對(duì)應(yīng)燃料包殼表面溫度最高點(diǎn)，這意味著堆芯水位降低會(huì)使燃料包殼表面?zhèn)鳠釔夯?，從而直接?dǎo)致PCT升高，因此及時(shí)注水提高堆芯水位能有效緩解事故。

根據(jù)圖5 可得，65 組工況中PCT 的歸一化值在一定區(qū)間內(nèi)變化，最大值為0.7504，最小值為0.7041，PCT的上限仍保持在1477K以下，表明堆芯有95%的置信度未發(fā)生損壞，且有較大的安全裕量，滿足安全驗(yàn)收準(zhǔn)則?？傮w來看，反應(yīng)堆堆芯在事故進(jìn)程中均能夠?qū)崿F(xiàn)有效冷卻，因此可以認(rèn)為在專設(shè)安全設(shè)施投入且滿足單一故障準(zhǔn)則的情況下，該船用壓水堆有能力在2000s 內(nèi)緩解48mm 中破口失水事故對(duì)堆芯的影響，實(shí)現(xiàn)堆芯的冷卻，確保反應(yīng)堆的安全。

3.2 敏感性分析

敏感性分析是在對(duì)數(shù)據(jù)結(jié)果開展不確定性分析的基礎(chǔ)上，通過度量各輸入?yún)?shù)對(duì)輸出的影響，識(shí)別出對(duì)目標(biāo)參數(shù)影響較大的關(guān)鍵參數(shù)。在統(tǒng)計(jì)學(xué)中，敏感性程度一般用相關(guān)性系數(shù)來衡量，常用的有Pearson和Spearman相關(guān)性系數(shù)。Spearman相關(guān)性系數(shù)是根據(jù)參數(shù)數(shù)值的位置排序來進(jìn)行相關(guān)性計(jì)算，與實(shí)際的數(shù)值大小無關(guān)，因此在輸入?yún)?shù)的數(shù)值量級(jí)相差較大或某個(gè)參數(shù)出現(xiàn)異常時(shí)仍然適用，且不用考慮采用歸一化等方法帶來的誤差。本文選用Spearman 相關(guān)性系數(shù)計(jì)算出各輸入?yún)?shù)與PCT 的相關(guān)性大小，從而篩選出MBLOCA 事故過程中對(duì)PCT影響較大的參數(shù)，該系數(shù)的表達(dá)式［12］為

式中，ρs為Spearman 相關(guān)性系數(shù)，Rxi為xi在輸入變量x中的大小排序，Ryi為yi在輸出變量y中的大小排序，n 為樣本數(shù)量。ρs取值范圍為-1～1，絕對(duì)值大小表示輸入?yún)?shù)與輸出參數(shù)的相關(guān)性強(qiáng)弱，絕對(duì)值越大表示相關(guān)性越強(qiáng)；符號(hào)表示相關(guān)性的正負(fù)，ρs＞0 表示兩個(gè)參數(shù)變化成正相關(guān)，即輸出參數(shù)隨輸入?yún)?shù)增大而增大，反之亦然，ρs=0表示兩個(gè)參數(shù)不相關(guān)。不確定輸入?yún)?shù)與PCT 的Spearman 相關(guān)性系數(shù)如圖6 所示（圖中編號(hào)對(duì)應(yīng)參數(shù)同表1）。

圖6 輸入?yún)?shù)與PCT之的Spearman系數(shù)

計(jì)算結(jié)果表明，對(duì)PCT 影響較大的輸入?yún)?shù)有：反應(yīng)堆運(yùn)行功率、低壓安注延遲時(shí)間、高壓安注延遲時(shí)間、破口面積、破口能量損失系數(shù)、低壓安注流量。其中，反應(yīng)堆運(yùn)行功率對(duì)目標(biāo)參數(shù)PCT的影響最大，這是因?yàn)樵谒プ円蜃硬蛔兊那闆r下，反應(yīng)堆運(yùn)行功率越大，堆芯衰變產(chǎn)生的熱量越多，PCT越高，即反應(yīng)堆運(yùn)行功率與PCT 變化成正相關(guān)。PCT 隨低壓安注延遲時(shí)間的增加而升高，這是因?yàn)榈蛪喊沧⒌难舆t投入會(huì)推遲堆芯再淹沒與冷卻的進(jìn)程，從而導(dǎo)致PCT升高。破口面積對(duì)應(yīng)的相關(guān)系數(shù)為負(fù)值，表明PCT 隨破口面積的增大而降低，實(shí)際上，由于低壓安注的流量遠(yuǎn)大于高壓安注流量，48mm中破口失水事故進(jìn)程中堆芯冷卻主要受低壓安注的影響，破口面積越大，系統(tǒng)壓力下降到低壓安注投入壓力整定值越快，低壓安注越早投入，從而使PCT越低。低壓安注流量越大，堆芯冷卻得更加充分，導(dǎo)致PCT越低。

4 事故緩解策略優(yōu)化

事故發(fā)生后，操縱員需要依據(jù)事故緩解策略及時(shí)正確的處置才能確保反應(yīng)堆的安全，因此，采用優(yōu)化的緩解策略可以增加反應(yīng)堆的安全裕量。根據(jù)上述分析，結(jié)合船用核動(dòng)力裝置的實(shí)際情況，可以采取提高低壓安注投入壓力和增加低壓安注流量?jī)煞N方式，提高反應(yīng)堆的事故緩解能力。選取6組優(yōu)化策略進(jìn)行對(duì)比分析，結(jié)果如表2所示，圖7為6 組策略對(duì)應(yīng)的燃料包殼表面溫度隨時(shí)間的變化曲線。

表2 優(yōu)化策略組合對(duì)比分析

比較表2 中1、2、3、4 組運(yùn)算結(jié)果可得，提高低壓安注投入壓力和增加低壓安注流量能夠有效降低PCT，增強(qiáng)反應(yīng)堆的事故緩解能力，但影響程度逐漸降低。比較1、5 組運(yùn)算結(jié)果，僅提高低壓安注投入壓力能夠顯著降低PCT；比較1、6 組運(yùn)算結(jié)果，僅增加低壓安注流量會(huì)導(dǎo)致PCT 略有降低；比較4、6 組運(yùn)算結(jié)果，將低壓安注流量增加相同程度，低壓安注投入壓力越高，PCT 越低；比較4、5 組運(yùn)算結(jié)果，將低壓安注投入壓力提高相同程度，低壓安注流量越大，PCT 越低，但影響遠(yuǎn)小于低壓安注投入壓力。綜上，考慮設(shè)計(jì)建造的成本及船用堆特性，該堆型可采取第3 組優(yōu)化策略來提高抵抗外界干擾的能力。

5 結(jié)語

針對(duì)船用核動(dòng)力裝置，采用最佳估算程序RELAP5/MOD3.2建模，并模擬冷管段中破口失水事故進(jìn)程，基于輸入不確定性傳播GRS方法開展不確定性量化計(jì)算，對(duì)結(jié)果進(jìn)行不確定性和敏感性分析，提出事故緩解優(yōu)化方案，得到以下主要結(jié)論：

1）GRS方法中采用Wilks公式的非參數(shù)統(tǒng)計(jì)方法進(jìn)行不確定性傳遞計(jì)算，有效地減少了抽樣運(yùn)算次數(shù)，基于Matlab程序開發(fā)的批量處理程序能有效提高分析效率。

2）在專設(shè)安全設(shè)施投入且滿足單一故障準(zhǔn)則的情況下，不確定性計(jì)算得到PCT最大歸一化值為0.7504，且堆芯在事故進(jìn)程中能夠?qū)崿F(xiàn)自我冷卻，表明堆芯有95%的置信度未發(fā)生損壞，滿足安全驗(yàn)收準(zhǔn)則

3）通過計(jì)算各輸入?yún)?shù)與PCT 之間的Spearman 相關(guān)性系數(shù)，得到對(duì)PCT 影響最大的參數(shù)為反應(yīng)堆運(yùn)行功率，即反應(yīng)堆事故前運(yùn)行功率越大，PCT 越高。低壓安注的延遲投入會(huì)推遲堆芯再淹沒與冷卻的進(jìn)程，且延遲時(shí)間長(zhǎng)短受人為因素的影響較大，操縱員正確判斷注水時(shí)機(jī)能有效避免燃料包殼表面?zhèn)鳠釔夯?，可為?yōu)化事故處置規(guī)程提供理論依據(jù)。

4）提高低壓安注投入壓力和增加低壓安注流量能夠增強(qiáng)反應(yīng)堆的事故緩解能力，但影響程度逐漸降低，考慮設(shè)計(jì)建造的成本及船用堆特性，通過對(duì)比分析選出合理可行的優(yōu)化策略。