基于抽樣的不確定性及敏感性分析的方法在核電廠水膜蒸發(fā)試驗誤差分析中的應用

扈本學，王　喆，王國棟，王章立，倪陳宵，張今朝，楊　萍（上海核工程研究設計院，上海　200233）

?

基于抽樣的不確定性及敏感性分析的方法在核電廠水膜蒸發(fā)試驗誤差分析中的應用

扈本學*，王喆，王國棟，王章立，倪陳宵，張今朝，楊萍
（上海核工程研究設計院，上海200233）

摘要：與傳統(tǒng)的誤差分析方法相比，基于抽樣的不確定性及敏感性分析具有較大的優(yōu)勢。本工作通過耦合DAKOTA程序和水膜蒸發(fā)試驗數(shù)據(jù)分析程序，開發(fā)了水膜熱態(tài)試驗誤差分析方法，計算得到了試驗目標參數(shù)水膜蒸發(fā)換熱乘子的不確定性范圍，并且分析了試驗測量參數(shù)的不確定性對蒸發(fā)換熱乘子不確定性的影響。計算結(jié)果表明，水膜入口流量、入口風速以及平板表面溫度是主要的不確定性來源。這為優(yōu)化試驗測量系統(tǒng)，減小試驗誤差提供了定量支持。該方法可以用于其他試驗誤差分析以及參數(shù)重要性分析。

關(guān)鍵詞：sobol方法；試驗誤差分析；敏感性分析；水膜蒸發(fā)試驗

試驗誤差分析是試驗數(shù)據(jù)分析中的重要一環(huán)，試驗的誤差范圍反映了試驗結(jié)果的精度，直接影響試驗的成敗與質(zhì)量。傳統(tǒng)的試驗誤差分析主要應用誤差傳遞的方法，根據(jù)直接參量、參數(shù)的誤差范圍，通過誤差傳遞公式推導，獲得試驗目標值的誤差。當試驗目標值的計算公式比較復雜時，采用該方法推導目標值的誤差較為困難。

不確定性分析和敏感性分析是分析復雜系統(tǒng)的重要工具［1］。不確定性分析研究輸入?yún)?shù)的不確定性通過模型傳播到模型輸出；而敏感性分析研究模型輸入的不確定性對模型輸出不確定性的貢獻率。主要分析思路：根據(jù)輸入?yún)?shù)的不確定性分布進行隨機抽樣，獲得輸入?yún)?shù)的不同組合（一個參數(shù)組合稱為一個抽樣工況），然后對每一個抽樣工況進行計算，獲得輸出參數(shù)，最后進行統(tǒng)計分析，獲得輸出參數(shù)的不確定性分布以及相關(guān)的敏感性度量參數(shù)。

該方法目前已經(jīng)在核電廠事故安全分析方面得到應用，主要用于驗收準則參數(shù)的不確定性分析及參數(shù)重要性評價。Gertman等人應用RELAP-3D程序，計算獲得了大破口失水事故峰值包殼溫度與各輸入?yún)?shù)之間的相關(guān)系數(shù)，定量評價了各個輸入?yún)?shù)的重要性［2-3］?；诓淮_定分析的方法，AREVA公司開展了大破口失水事故后核電廠安全殼壓力響應的不確定性分析，并利用計算得到的相關(guān)系數(shù)，分析了影響安全殼壓力的各輸入?yún)?shù)的重要程度［4］。近年來，在反應堆工程領(lǐng)域，我國業(yè)內(nèi)人士也陸續(xù)開展了基于抽樣的不確定性及敏感性分析的工作，包括：熔鹽堆參數(shù)的不確定性分析［5］、大破口失水事故后質(zhì)能釋放參數(shù)敏感性分析［6］、設計基準事故后安全殼壓力響應的敏感性分析［7］以及鈉冷快堆事故不確定性分析等［8］。

本工作將基于隨機抽樣不確定性分析及敏感性分析的方法引入試驗的誤差分析中，以水膜蒸發(fā)換熱乘子為試驗目標值，通過對試驗的直接測量參數(shù)進行抽樣，計算獲得試驗目標值，然后進行統(tǒng)計分析，獲得試驗目標值的誤差帶，并利用計算的Sobol敏感度系數(shù)來判斷試驗目標值的主要誤差來源。

1　基本理論

1.1容許區(qū)間與Wilks公式

應用隨機抽樣方法進行不確定性分析時，往往需要確定抽樣次數(shù)，使得輸出參數(shù)滿足特定置信度和特定概率的要求。

假設輸出參數(shù)為一系列輸入?yún)?shù)的函數(shù)：y=f（x1，x2，…，xN），其概率分布函數(shù)為g（y）。進行N次抽樣，可以計算得到因變量y的一個樣本｛y1，y2，…，yN｝。則可以確定兩個函數(shù)L=L（y1，y2，…，yN），U=U（y1，y2，…，yN），使得：

其中，β為置信水平，γ為概率比率，［L，U］為容許區(qū)間。公式（1）表示的意義如下：隨機變量y的m次計算結(jié)果有一定分布，即在置信度為β時，計算結(jié)果落入?yún)^(qū)間［L，U］的次數(shù)占總計算次數(shù)m的概率比率大于γ，對于安全分析，希望β和γ的值盡可能高。對于給定的β和γ，確定隨機抽樣的次數(shù)m后，便可以確定滿足式（1）的容許區(qū)間［L，U］。

對于函數(shù)y=f（x1，x2，…，xN），用Wilks公式［9］，即采用隨機抽樣加排序統(tǒng)計方法來確定相應的容許區(qū)間。假設y1，y2，…，yN為N次獨立的隨機函數(shù)y的輸出，其概率密度分布函數(shù)g（y）是連續(xù)的。將y1，y2，…，yN按升序排列，定義y（k）為序列中的第k個值，可以得出：y（1）= min（yk），y（N）=max（yk）。

若取最大值和最小值作為雙側(cè)容許區(qū)間，其置信水平為：

若取最大值作為單側(cè)容許區(qū)間，則其置信水平為

對于自變量與因變量的任意函數(shù)關(guān)系，Wilks公式都可以保證容許區(qū)間滿足相應的β 和γ。

根據(jù)Wilks公式，根據(jù)輸入?yún)?shù)的概率密度分布，通過隨機抽樣可以確定93組狀態(tài)點，通過函數(shù)關(guān)系式可以計算得到輸出函數(shù)y=f（x1，x2，…，xN）的樣本（y1，y2，…，y93）。假設樣本的最大值為ymax，最小值為ymin，則滿足兩個95%的雙側(cè)容許區(qū)間為［ymin，ymax］。其物理意義為在置信度為95%時，輸出函數(shù)y所有計算結(jié)果中落入?yún)^(qū)間［ymin，ymax］的概率大于95%。因此，對于水膜熱態(tài)試驗誤差分析，應用93次隨機抽樣，可以獲得滿足兩個95%的試驗目標參數(shù)的上下限。

1.2Sobol方法

Sobol［10，11］方法是最代表性的全局參數(shù)敏感性分析方法，基于模型分解思想，得到輸入?yún)?shù)一、二次甚至更高次的敏感度。國內(nèi)學者已開始將Sobol方法用于參數(shù)的重要性分析。文獻［12］通過比較各因素對換熱網(wǎng)絡的一階靈敏度及總靈敏度，甄別出對換熱網(wǎng)絡系統(tǒng)運行影響較大的因素以及影響較小但存在交互影響的因素。文獻［13］用Sobol方法開展了水文模型參數(shù)敏感性分析。

Sobol方法一共包含兩個指數(shù)，主要敏感度系數(shù)Si和總敏感度系數(shù)Ti。其主要理論依據(jù)如下［14］：

假定y=f（x）（x1，x2，…，xn），并且xi服從［0，1］均勻分布，且f2（x）可積，則函數(shù)f（x）（可以分解為：

則模型的方差D也可以分解為單個參數(shù)和每個參數(shù)的影響：

對上式歸一化，令：

稱為Sobol敏感度系數(shù)。Si稱為一階敏感度系數(shù)，又稱為主要敏感度系數(shù)；S1，2，…，n稱為n階敏感度系數(shù)?？偯舾卸认禂?shù)的計算公式為：

其中，主要敏感度系數(shù)Si體現(xiàn)了單個輸入?yún)?shù)的不確定性對輸出參數(shù)方差的貢獻程度，全部敏感度系數(shù)Si體現(xiàn)了單個輸入?yún)?shù)不確定度以及該參數(shù)與其他參數(shù)的相互作用對輸出參數(shù)方差的貢獻程度。當Si與Ti的值差距較大時，說明因素的交互作用明顯。

2　水膜熱態(tài)試驗

2.1試驗簡介

由于CAP1400安全殼尺寸較AP1000有所增大，反應堆核功率和設計基準事故（Design Basis Accident，簡稱DBA）后安全殼的峰值壓力也有所增大，為了在更大范圍內(nèi)充分驗證CAP1400安全殼分析程序水膜蒸發(fā)關(guān)系式的適用性，設立了重大專項非能動安全殼冷卻系統(tǒng)（Passive Containment Cooling System，簡稱PCCS）水膜熱態(tài)試驗。非能動安全殼冷卻系統(tǒng)水膜熱態(tài)試驗的核心目標是在CAP1400非能動安全殼冷卻系統(tǒng)系統(tǒng)運行參數(shù)范圍內(nèi)驗證目標經(jīng)驗關(guān)系式，并且通過試驗驗證水膜蒸發(fā)傳熱傳質(zhì)關(guān)系式包絡因子的保守性。

圖1　水膜熱態(tài)試驗示意圖Fig.1 Sketch of heated flat plate test

試驗臺架如圖1所示。試驗本體由加熱平板、兩個側(cè)壁以及透明的有機玻璃罩組成。加熱平板采用不銹鋼平板，并安裝玻璃罩形成矩形風道和水道。矩形通道的寬度為1.183 m，高度為0.285 m，總長度為7.3 m。板后鋪設了加熱銅管，通入流動的導熱油為試驗提供所需熱源，被加熱銅管覆蓋的平板長度為5 m。試驗本體的表面粘有熱流密度計，用以測量平板表面溫度以及熱流密度。試驗中需要測量的參數(shù)有：入口冷卻水和空氣的流量及溫度，不同位置處水膜的流速、溫度、寬度（覆蓋率）、厚度，不同位置處空氣的流速、溫度、濕度，不同位置處的平板溫度、表面熱流密度，平板傾斜角度等。

2.2目標關(guān)系式驗證

平板水膜蒸發(fā)試驗的目的是驗證水膜與空氣之間的傳質(zhì)關(guān)系式符合相關(guān)目標經(jīng)驗關(guān)系式。對于目前應用的WGOTHIC分析程序，應用以下的方法來預測水膜蒸發(fā)傳質(zhì)系數(shù)。

當主流氣體中的蒸汽濃度和液體表面的蒸汽濃度不相等時，它們之間就會發(fā)生質(zhì)量傳遞。蒸汽的濃度梯度指的是主流氣體中和液體表面的蒸汽分壓的差值。如果主流氣體中的蒸汽濃度高于液體表面的蒸汽濃度，就會發(fā)生冷凝；如果液體表面的蒸汽濃度高于主流氣體中的蒸汽濃度，就會發(fā)生蒸發(fā)。

類似于對流換熱過程，可以引出以下公式：

式中，Gevap.是蒸發(fā)/冷凝質(zhì)量流量，kg是傳質(zhì)系數(shù)，Mv是蒸汽的分子量，Pv，surf.是液體表面的蒸汽分壓，Pv，bulk是主流氣體中的蒸汽分壓。

舍伍德數(shù)（Sh）是反映傳質(zhì)系數(shù)的無量綱數(shù)，類似對流換熱的努賽爾數(shù)（Nu）定義為：

式中，R是摩爾氣體常數(shù)；T是邊界層溫度，取主流體溫度和液膜溫度的平均值；L是特征長度；P是總壓；Plm，air是空氣對數(shù)平均分壓；Dv是蒸汽在空氣中的擴散系數(shù)。Plm，air、Dv分別由式（11）和（12）計算：

式中，Pair，bulk是主流氣體中空氣分壓，Pair，surf.是液膜表面處空氣分壓。

式中，T（°F）為主流溫度和液膜溫度的平均值，P（psia）為主流氣體的總壓。對于傳熱關(guān)系式，用Sh代替Nu、用施密特數(shù)（Sc）代替普朗特數(shù)（Pr）可以轉(zhuǎn)換成傳質(zhì)關(guān)系式，得到以下關(guān)系式：

式中，Sc=v/Dv，v為運動粘性系數(shù)；Nu數(shù)采用以下公式計算：

逆浮升力方向?qū)α鲹Q熱關(guān)系式：

順浮升力方向?qū)α鲹Q熱關(guān)系式：

式中，

式中，Re為雷諾數(shù)，Gr為表征自然對流強度的格拉曉夫數(shù)。

對于本試驗，浮升力方向與空氣運動方向相同，屬于順浮升力方向混合對流換熱，因此，選取順浮升力方向?qū)α鲹Q熱關(guān)系式。

應用試驗測量的水膜蒸發(fā)率，并應用公式（8）～公式（12），可以得到試驗的傳質(zhì)舍伍德數(shù)Shmeas。而應用試驗的Sc、Pr、Re數(shù)值，應用公式（13）可以獲得預測得到一個舍伍德數(shù)Shpred。

定義蒸發(fā)因子為：Shmeas/Shpred。對于本試驗，核心的目標是獲得每個試驗工況的蒸發(fā)因子。進行誤差分析即獲得蒸發(fā)因子的不確定性范圍，因此進行不確定性分析和敏感性分析時目標值（輸出參數(shù)）為蒸發(fā)因子。

3　分析方法

DAKOTA程序是美國圣迪亞國家實驗室開發(fā)的最優(yōu)化、參數(shù)估計、不確定性分析及敏感性分析程序。該程序擁有強大的接口能力，可與其它計算程序耦合連接［12］。對于水膜熱態(tài)試驗的誤差分析，應用DAKOTA程序耦合水膜熱試驗數(shù)據(jù)分析程序?qū)崿F(xiàn)蒸發(fā)因子的不確定性及敏感性分析。

應用DAKOTA程序耦合水膜蒸發(fā)因子計算程序，開發(fā)了水膜熱態(tài)試驗不確定性分析及敏感性分析的工具。如圖2所示：DAKOTA程序根據(jù)單個試驗工況的試驗測量參數(shù)（抽樣參數(shù)）及其不確定性分布進行隨機抽樣，獲得了N個輸入?yún)?shù)的組合，每個參數(shù)組合代表一個計算工況；然后DAKOTA程序調(diào)用水膜熱態(tài)試驗數(shù)據(jù)程序分別計算獲得每一個抽樣工況的蒸發(fā)因子，并將其返回DAKOTA程序；最終應用DAKOTA程序統(tǒng)計分析計算結(jié)果。

應用抽樣方法進行敏感性分析及不確定分析時，首先需要確定抽樣變量及其不確定性分析，根據(jù)試驗測量結(jié)果，一共選取13個測量參數(shù)作為抽樣參數(shù)，其測量誤差見表1。表1中最后一列給出了誤差類型，“絕對”表示給出是參數(shù)的實際誤差，而“相對”表示給出的是參數(shù)的實際誤差/名義值。通過抽樣獲得各個參數(shù)誤差，然后根據(jù)誤差類型，將該值或該值與名義值的乘積與名義值相加即可得到該抽樣工況的參數(shù)值。

根據(jù)Wilks公式，對于不確定性分析，對每個試驗工況應用蒙特卡羅抽樣93次，計算結(jié)果獲取的蒸發(fā)因子的最大值和最小值就是獲得滿足兩個95%的試驗目標參數(shù)的上下限值。而對于蒸發(fā)因子的敏感性分析，根據(jù)要求［15］，抽樣次數(shù)應為N×（M+2），N至少為幾百乃至上千，M為抽樣變量的數(shù)目。分析中應用蒙特卡羅抽樣，抽樣次數(shù)為5 000×15，統(tǒng)計獲得輸入?yún)?shù)（試驗測量參數(shù)）的總的Sobol敏感度系數(shù)Ti。

圖2　程序耦合流程Fig.2 Flow Chart of code coupling

表1　測量參數(shù)的測量不確定性Table 1 Uncertainty of measured parameters

4　計算結(jié)果

圖3給出了水膜熱態(tài)試驗5個試驗工況的蒸發(fā)因子的不確定性分析結(jié)果。圖中最大值、最小值對應該工況下93個抽樣計算蒸發(fā)因子的最大值和最小值，而名義值代表的是應用試驗各測量參數(shù)的名義值進行計算獲得的蒸發(fā)因子?？梢缘贸觯?/p>

（1）各個試驗工況蒸發(fā)因子的名義值偏離1.0很小，最大偏差大約為0.2，說明目標關(guān)系式可以較好的模擬水膜蒸發(fā)換熱。

（2）考慮各工況測量參數(shù)的不確定性以后獲得的蒸發(fā)因子的上下限值與名義值的最大相對偏差低于10%，說明試驗目標值的不確定性較小，試驗精度較高。

圖3　蒸發(fā)因子的不確定性上下限值Fig.3 Upper and lower limits of evaporation ratio uncertainty

圖4～圖8分別給出了5個試驗工況中各測量參數(shù)的總的Sobol敏感度系數(shù)Ti的值，可以得出：雖然影響各個工況的Ti值不盡相同，但綜合評價Ti較大的參數(shù)為：水膜進、出口流量、5個平板溫度測點以及入口風速。工況2中入口風速的不確定性對蒸發(fā)因子影響很小的原因是該工況風速的名義值大約為1.9 m·s-1，位于自然對流區(qū)，風速對程序預測的蒸發(fā)量沒有影響，因此對蒸發(fā)因子的值影響很小。

圖9給出了5個工況各測量參數(shù)的總的Sobol敏感度系數(shù)Ti的平均值。根據(jù)圖8的結(jié)果，評價水膜熱態(tài)試驗誤差中最重要的來源（超過1%）是：水膜入口流量、入口風速以及平板溫度測點1～5。

圖4　試驗工況1的總敏感度系數(shù)Fig.4 Total sensitivity coefficients of Test 1

圖5　試驗工況2的總敏感度系數(shù)Fig.5 Total sensitivity coefficients of Test 2

圖6　試驗工況3的總敏感度系數(shù)Fig.6 Total sensitivity coefficients of Test 3

圖7　試驗工況4的總敏感度系數(shù)Fig.7 Total sensitivity coefficients of Test 4

圖8　試驗工況5的總敏感度系數(shù)Fig.8 Total sensitivity coefficients of Test 5

圖9　5個試驗工況的總敏感度系數(shù)平均值Fig.9 Average total sensitivity coefficients of 5 tests

5　結(jié)論

本工作應用DAKOTA程序耦合水膜熱態(tài)試驗處理程序，開展了基于抽樣的水膜熱態(tài)試驗蒸發(fā)因子的不確定性及敏感性統(tǒng)計分析。分析結(jié)果表明：水膜熱態(tài)試驗的誤差控制較好，水膜熱態(tài)試驗獲得蒸發(fā)換熱因子的不確定性較小，試驗精度較高；水膜熱態(tài)試驗的誤差主要來源為水膜入口流量、入口風速以及平板溫度測量。該方法提供了一種用于確定試驗數(shù)據(jù)的不確定性及不確定性來源的方法，分析結(jié)果可以為試驗測量系統(tǒng)的優(yōu)化提供支撐。

參考文獻

［1］J. C. Helton and F. J. Davis. Illustration of Sampling-Based Methods for Uncertainty and Sensitivity Analysis［J］. Risk Analysis，2002，22（3）：591-622.

［2］Idaho National Laboratory. Uncertainty Analysis for RELAP5-3D［R］. Idaho：Idaho National Laboratory，2011.

［3］Idaho National Laboratory. Uncertainty Analysis of RELAP5-3D［R］. Idaho：Idaho National Laboratory，2012.

［4］Jamal M. Abdelghany，Robert P. Martin. Uncertainty Analysis for Containment Response of U.S. EPR？Reactor to Large Break Loss-of-Coolant Accidents［C］. Sanhigeo：Proceedings of ICAPP '10，2010：1395-1409.

［5］王成龍，Lin-wen Hu，秋穗正，等.可移動式熔鹽冷卻高溫堆（TFHR）熱工水力特性不確定性分析［C］.北京：第十四屆熱工流體會議論文集，2015.

［6］扈本學，王喆，王國棟，等.應用DAKOTA程序耦合WCOBRA/TRAC程序進行大破口失水事故質(zhì)能釋放參數(shù)敏感性分析［C］.北京：第十四屆熱工流體會議論文集，2015.

［7］王國棟，王喆，扈本學，等.應用DAKOTA程序耦合WGOTHIC程序進行安全殼壓力響應敏感性分析［C］.北京：第十四屆熱工流體會議論文集，2015.

［8］岳倪娜，馬在勇，蔡容，等.鈉冷快堆事故不確定性分析［C］.北京：第十四屆熱工流體會議論文集，2015.

［9］Wilks S S. Determination of sample sizes for setting tolerance limits［J］. The Annals of Mathematical Statistics，1941，12 （1）：91-96.

［10］Sobol IM. Sensitivity analysis for non-linear mathematical models［J］. Math Modeling Comput Exp，1993（1）：407-414.

［11］Homma T，Saltelli A. Importance measures in global sensitivity analysis of nonlinear models［J］. Reliability Engineering and System Safety，1996，52（1）：1-17.

［12］張紅晶，吉海濤，譚世語，等.基于Sobol'法的換熱網(wǎng)絡全局靈敏度分析［J］.世界科技研究與發(fā)展，2012，34（6）：916-919.

［13］任啟偉，陳洋波，周浩瀾，等.基于Sobol法的TOPMODEL模型全局敏感性分析［J］.人民長江，2010，41（19）：91-107.

［14］V. Gregory Weirs，James R. Kamm，Laura P. Swiler，et al. Sensitivity analysis techniques applied to a system of hyperbolic conservation laws［J］. Reliability Engineering and System Safety，2012，107（12）：157-170.

［15］Dakota，a multilevel parallel object-oriented framework for design optimization，parameter estimation，uncertainty quantification，and sensitivity analysis，akota user's manual［M］. Sandia National Laboratories，Version 5.3.1，2013.

Aplication of the Sampling-Based Uncertainty and Sensitivity
Analysis Method in Error Analysis of Heated Flat Plate Test at NPP

HUBenxue，WANGZhe，WANGGuodong，WANGZhangli，NIChenxiao，ZHANGJinzhao，YANGPing
（Shanghai Nuclear Researchand Design Institute，Shanghai200233，China）

Abstract:Compared with traditional error analysis methods，uncertainty and sensitivity analysis based on sampling has great advantages. Method for error analysis of heated flat plate test has been developed by coupling DAKOTA code and data evaluation program. Uncertainty of evaporation ratio has been calculated and effects of measured parameters uncertainty on evaporation ratio have been analyzed. The calculated results show that inlet film flow rate，inlet air flow velocity and plate surface temperatures are the main sources of the uncertainty，which provide support for measure system optimization to reduce the error range. This method could be applied to the error analysis for other tests and parameters importance analysis.

Keywords:Sobol method；test error analysis；sensitivity analysis；heated flat plate test

中圖分類號：TL364

文章標志碼：A

文章編號：1672-5360（2016）01-0084-06

收稿日期：2015-10-08修回日期：2015-11-23

基金項目：國家科技重大專項，項目編號2011ZX06002-0056

作者簡介：扈本學（1986—），男，山東濰坊人，工程師/博士，反應堆熱工，現(xiàn)主要從事反應堆安全分析及非能動安全殼分析工作

*通訊作者：扈本學，E-mail：hubenxue@snerdi.com.cn