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

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

?

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

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

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

關鍵詞：sobol方法；試驗誤差分析；敏感性分析；水膜蒸發試驗

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

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

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

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

1　基本理論

1.1容許區間與Wilks公式

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

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

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

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

若取最大值和最小值作為雙側容許區間，其置信水平為：

若取最大值作為單側容許區間，則其置信水平為

對于自變量與因變量的任意函數關系，Wilks公式都可以保證容許區間滿足相應的β 和γ。

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

1.2Sobol方法

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

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

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

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

對上式歸一化，令：

稱為Sobol敏感度系數。Si稱為一階敏感度系數，又稱為主要敏感度系數；S1，2，…，n稱為n階敏感度系數。總敏感度系數的計算公式為：

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

2　水膜熱態試驗

2.1試驗簡介

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

圖1　水膜熱態試驗示意圖Fig.1 Sketch of heated flat plate test

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

2.2目標關系式驗證

平板水膜蒸發試驗的目的是驗證水膜與空氣之間的傳質關系式符合相關目標經驗關系式。對于目前應用的WGOTHIC分析程序，應用以下的方法來預測水膜蒸發傳質系數。

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

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

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

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

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

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

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

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

逆浮升力方向對流換熱關系式：

順浮升力方向對流換熱關系式：

式中，

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

對于本試驗，浮升力方向與空氣運動方向相同，屬于順浮升力方向混合對流換熱，因此，選取順浮升力方向對流換熱關系式。

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

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

3　分析方法

DAKOTA程序是美國圣迪亞國家實驗室開發的最優化、參數估計、不確定性分析及敏感性分析程序。該程序擁有強大的接口能力，可與其它計算程序耦合連接［12］。對于水膜熱態試驗的誤差分析，應用DAKOTA程序耦合水膜熱試驗數據分析程序實現蒸發因子的不確定性及敏感性分析。

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

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

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

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

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

4　計算結果

圖3給出了水膜熱態試驗5個試驗工況的蒸發因子的不確定性分析結果。圖中最大值、最小值對應該工況下93個抽樣計算蒸發因子的最大值和最小值，而名義值代表的是應用試驗各測量參數的名義值進行計算獲得的蒸發因子。可以得出：

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

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

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

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

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

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

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

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

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

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

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

5　結論

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

參考文獻

［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程序進行大破口失水事故質能釋放參數敏感性分析［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'法的換熱網絡全局靈敏度分析［J］.世界科技研究與發展，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—），男，山東濰坊人，工程師/博士，反應堆熱工，現主要從事反應堆安全分析及非能動安全殼分析工作

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