校準條件下的數值模擬不確定度量化方法

馬智博，殷建偉，李海杰，劉 全

（1.北京應用物理與計算數學研究所，北京 100088；2.長沙機電產品研究開發中心，湖南長沙 410100）

文章編號：1001?246X（2015）05?0514?09

校準條件下的數值模擬不確定度量化方法

馬智博1，殷建偉1，李海杰2，劉 全1

（1.北京應用物理與計算數學研究所，北京 100088；2.長沙機電產品研究開發中心，湖南長沙 410100）

當存在眾多不確定輸入因素時，不確定度的傳遞分析往往導致對數值模擬不確定度的過高估計.利用校準行為能夠消減系統級數值模擬中認知不確定度的客觀機制，提出一個綜合利用已有系統級試驗對比信息和新增建模與模擬傳遞信息的不確定度量化方法，結合一個虛擬試驗的例子對該方法進行展示和驗證.

數值模擬；不確定度量化；校準；可靠性認證

0 引言

在投入應用之前，建模與模擬（Modeling＆Simulation，M＆S）技術需要經過驗證（Verification）、確認（Validation）和認可（Accreditation）過程，以評估其可信度并獲取工程應用所需的資質［1－2］.事實上，即便經過充分的驗證和確認，也難以完全發現并徹底解決M＆S中的所有問題，再加上不可避免的離散誤差，對復雜物理過程的數值模擬往往會產生相對于真值的系統性偏離，因而常常通過校準（Calibration）措施給予糾偏以提高模擬結果與試驗結果的一致性.用經過校準的模擬程序對新的復雜對象進行模擬計算，是現實中普遍存在并將長期存在的工程預測方法.當不能進行新對象的系統級試驗時，對校準后數值模擬的預測結果進行不確定度量化，是可靠性認證與評估所面臨的一個重要科學問題［3］.

認證的對象為抽象的事物類（如型號），評估的對象為具體的事物例（如產品）［4］.型號經過認證后獲得定型，定型產品經過庫存和服役等過程后，可能發生狀態上的改變，往往需要分期評估.為全面考察各種工程因素的影響，新型號認證一般會安排較多的系統級試驗，試驗數據可用于數值模擬的系統級校準、確認和不確定度量化.當不能進行新對象的系統級試驗時，需要針對改變后的設計參數和新增加的工程因素進行數值模擬及其模擬結果的不確定度量化.

如果不能進行新對象的系統級試驗，數值模擬不確定度量化的信息主要來自兩個方面：其一是對比信息，即通過原型對象的系統級試驗數據與數值模擬的對比而得到的不確定度信息，利用模型空間中不確定度從確認域到應用域的外推技術，可以將該信息轉換為能夠直接反映新對象系統級數值模擬的不確定度信息［4］；其二是傳遞信息，即借助總體計算程序和擾動分析方法，將較低層級建模與模擬不確定度信息傳遞到新對象系統級數值模擬的不確定度信息［5］.

可靠性分析往往受困于信息缺乏，多源信息的整合因此變得迫切.雖然對比信息和傳遞信息分別來自不同的認知途徑，但反映的是同一個問題，將二者進行整合，既有現實需求，又有科學上的合理性和可行性.

基于單信息源的不確定度量化已有較好的方法基礎.Oberkampf根據數值模擬結果以及試驗結果的均值和標準差，給出了概率框架下基于對比分析的數值模擬誤差量化方法［1］，Helton總結了認知不確定度傳遞分析的敏感度方法和蒙特卡羅抽樣方法［6］，劉全用非嵌入多項式混沌方法進行了炸藥狀態方程模型參數不確定度對爆轟計算結果不確定度的傳遞分析［7］.將對比信息從確認域外推到應用域的不確定度量化方法也被提出并逐漸投入使用［4，8］，而有關這兩類信息的整合，仍是目前面臨的技術瓶頸.

數值模擬的不確定度主要為認知不確定度，宜用區間理論來表達并進行信息整合，不確定度量化應當遵循真值覆蓋原則和最小化原則［4］.如果量化結果僅利用外推后的對比信息，則會因為沒有充分考慮新對象建模和模擬所新增加的認知缺陷，容易導致不確定度估計結果偏小而違背真值覆蓋原則，最終增加不合格對象被接收的風險，如果把外推后的對比信息和不確定度的整體傳遞信息直接相加或僅僅利用整體傳遞信息，均可能使不確定度估計結果被不合理地放大而違背最小化原則，從而增加合格對象被拒收的風險.

根據認知科學中對比信息和傳遞信息之間的容斥關系以及不確定度量化的基本原則，本文將經過外推的系統級對比信息作為不確定度的基本組成部分，將新對象建模和模擬所新增加的不確定度傳遞到系統級的計算結果中，作為不確定度的增量部分，然后將二者進行相加作為新對象數值模擬的不確定度，并結合一個黎曼問題及其虛擬試驗對該不確定度量化方法進行展示和驗證.

1 數值模擬校準及其技術特征

提高計算結果與試驗結果一致性的途徑有兩個，其一是全面提高認知能力，降低數值模擬中每個技術環節的認知不確定度，這也是數值模擬發展的理想途徑，其二是基于現有認知能力，人為地讓不同技術環節所產生的認知誤差相互補償.數值模擬校準在較大程度上依賴于誤差補償機制，但這種補償作用只在較小的模型范圍內有效，隨著模擬對象遠離那些用以校準的試驗對象，數值模擬不確定度往往增加很快.

本文提及的模型包括實體模型和物理模型，前者指具體的特殊工程對象，包括部件的材料選用、密度與質量、形狀與尺寸、動作過程的初邊條件等工程設計參數，后者指抽象的一般物理規律，包括描述材料力學屬性的本構關系和狀態方程、描述相關物理過程（如爆轟、湍流和傳熱等）的模型方程以及描述物質運動普遍規律的三大守恒方程等.實體建模的不確定度往往決定于靜態測量的不確定度，物理建模可能涉及更多的動態測量以及基于測量結果和近似理論的逆向運算，因此，物理建模所引起的數值模擬不確定度一般較大.

結合試驗進行的校準行為，主要體現在對模型形式、計算方法和選定的校準參數進行調整并把它們固定下來［2］，校準參數一般來自計算參數、物理參數或可調參數.可調參數可理解為缺乏明確數學和物理含義或缺乏確切賦值信息但對最終數值模擬結果有明顯影響的參數.通過充分的驗證和確認，可有效減少可調參數的數量，提高數值模擬的預測能力，但由于數值計算對應非零的時間和空間離散尺度，離散誤差不可避免，再加上模型形式的近似性和計算方法在模擬復雜問題時能力的欠缺，可調參數的存在往往難以避免［9］.

校準過程一般基于某一范圍內的實體模型，本文將實體模型空間中的這個范圍稱為校準域.校準之后所進行的確認活動，也有對應的確認域，其中模擬結果和試驗結果的對比分析，被用于確認域內數值模擬的不確定度量化.

校準過程完成以后，計算程序及其相關的模型參數即被固化，但固化通常具有相對性和階段性.相對性主要表現在描述物理規律的模型參數在取值上與特定計算方法和計算參數有關，并非完全基于對物理的認識，而階段性表現在固化版本通常是暫時的，隨著數值模擬技術的發展，模型和參數也會隨之改變.

演繹和歸納是認識事物的兩個基本方法，前者基于特定條件和普遍規律來推斷事物的發展結果，是不確定度傳遞分析的基礎，后者通過實踐和觀察來總結事物的發展結果，是不確定度對比分析的基礎.對比信息直接來自系統層級的實踐活動，因此可信度更高，在進行信息整合時，對比信息的權重因子一般要遠大于傳遞信息.

我們認為，經歷校準后的數值模擬具有如下特點：

1）校準域內系統級數值模擬的不確定度能夠得到有效消減，但隨著實體模型遠離校準域，誤差補償作用逐漸減弱，新對象模擬結果對真值的背離會逐漸增大；

2）確認域內由對比分析得到的不確定度可外推至應用域，但外推結果不能充分反映新對象建模和模擬所導致的新增不確定度；

3）系統級以下物理建模和數值計算的不確定度并不因為校準而得到必然的消減，由傳遞分析得到的不確定度往往較大；

4）如果對比信息和傳遞信息均反映同一數值模擬的不確定度，進行信息整合時，權重分配應以對比信息為主；

5）當確認域外的新對象不能進行試驗時，數值模擬不確定度存在兩個獨立的信息來源：其一是系統級確認域內來自對比的不確定度信息，其二是系統級以下新增的建模和模擬不確定度所傳遞到系統級的不確定度信息，二者具有可加性，可基于區間理論進行信息整合.

2 校準條件下的不確定度量化

需要解決的關鍵問題之一，是如何整合對比信息和傳遞信息，并使量化方法符合真值覆蓋原則和最小化原則［4］.

實物試驗可消減認知不確定度，但不能消減偶然不確定度.對偶然不確定度，可根據各自的信息量來分配信息整合的權重，然后用概率框架下的貝葉斯理論或其它技術途徑進行對比信息和傳遞信息的整合［10－11］.數值模擬主要產生認知不確定度，只要對比信息存在，一般可把傳遞信息的權重賦零，主要使用對比信息進行不確定度量化.

當沒有新對象的系統級試驗時，無直接的對比信息可利用，不確定度量化的可能方案有：

1）僅僅利用傳遞信息；

2）僅僅利用外推后的對比信息；

3）同時利用以上兩種信息.

第一方案未利用原有的系統級對比信息，由傳遞給出的認知不確定度不能得到有效消減，傳遞過程產生的不確定度難以量化，低層級數值計算造成的不確定度也難以傳遞到系統層級.第二方案可能漏掉新對象建模和模擬所新產生的不確定度.第三方案有更加合理的思想基礎，但需要設計一套恰當的實施方法，最大限度地減少上述問題造成的影響.

基于第三方案，本文建議，首先將新對象由傳遞得到的數

傳遞信息一般由三部分組成，結合式（1），有如下近似表達式

其中1Upropagation，2Upropagation，3Upropagation分別表示實體建模、物理建模和數值計算三個技術環節所傳遞到系統級的數值模擬不確定度.

新對象數

綜上所述，可據如下步驟進行數值模擬不確定度的量化：

1）結合校準試驗，進行數值模擬校準，然后完成對計算程序和相關參數的固化；

2）結合確認試驗，進行系統級確認域內的對比分析和不確定度量化；

3）根據不確定度隨實體模型參數的變化規律，用確認域內的不確定度信息外推新對象數值模擬的不確定度；

5）根據公式（3）得到新對象數值模擬的總體不確定度.

3 基于對比分析的不確定度量化

試驗對象的加工隨機性和測試隨機性都會使試驗結果產生偶然不確定度，但該不確定度并非由數值模擬產生，因此，基于對比分析的不確定度量化方法，要恰當地構造統計量，屏除這類不確定度的影響，同時有效地揭示真正的數值模擬不確定度信息.

在工程實踐中，基于試驗的對比分析有兩種情形：

1）一對一，即一次數值模擬對應一次試驗；

2）一對多，即一次數值模擬對應多次試驗.

在一對一情形中，需要對每次試驗的對象進行單獨測量，然后將測量結果用于實體建模和物理建模，每個試驗對象均有相應的數值模擬結果和試驗結果，在扣除測試誤差后，兩者之間的差別能夠直接地反映試驗對象的數值模擬誤差.如果試驗和模擬次數均為n，可按下式提取數值模擬的不確定度：

一對多情形往往對應于重復性試驗，僅基于設計方案的公稱參數或實際加工的中值參數進行實體建模和物理建模，最終只產生一套數值模擬結果.雖然每一個試驗對象都按照同樣的設計方案進行加工制造，但產品加工和試驗測試的隨機性會導致試驗結果的隨機性.為了屏除隨機性干擾信息，本文建議從數值模擬結果與試驗均值之間的差異中挖掘數值模擬的不確定度信息，具體公式

本文取置信度β＝0.95，t（1－β）/2，ν為自由度為ν的t分布的（1－β）/2分位數，ν＝n－1．

當所有試驗對象均處于某設計點附近且n足夠大時，式（4）、（7）給出的結果趨于一致.

4 實例演示

4.1 問題描述

圖1所示的實體模型對應一維黎曼問題，共有兩個物質區域，左端區域相當于高壓氣體，其初始時刻的壓強、密度、內能和速度分別為p1，ρ1，e1，v1，右端區域相當于重金屬材料，其初始時刻的壓強、密度、內能和速度分別為p2，ρ2，e2，v2．兩區在零時刻的初速度均為零，從零時刻開始，在物質界面處產生一個向左傳播的膨脹波和向右傳播的沖擊波，金屬材料在高壓氣體的推動下向右運動，本例所關心的系統級輸出量是10μs時物質界面移動的距離D.

圖1 實體模型示意圖Fig.1 Sketchmap of entitymodel

物理模型主要涉及兩種材料的狀態方程，左端的氣體選用理想氣體狀態方程p＝（γ1－1）ρe，音速方程c右端的金屬也被視作流體，選用凝聚介質實用狀態方程p＝（ρ－ρ0）＋（γ2－1）ρe，音速方程 ，其中γ1和γ2分別為兩種材料的多方指數，c0為初始音速.

本例將上述實體模型看作一個產品型號，并假定庫存老化僅導致兩個物理參數發生改變，其一是氣體的初始內能e1（初始壓強與初始內能相關，不作為獨立的模型變量），其二是金屬的多方指數γ2．因二者均為庫存時間t的函數，可基于庫存時間建立實體模型的模型空間，并在該模型空間中定義校準域、確認域和應用域.

要解決的問題是，當不能對庫存產品進行新的系統級試驗時，如何基于已有的0年～50年庫存期產品的系統級試驗結果及其對應的數值模擬結果，量化80年庫存期產品數值模擬預測的不確定度.

4.2 試驗設計

模型校準、模型確認以及對不確定度量化方法的檢驗將涉及試驗，該例存在系統級輸出量的精確解D?，因此擬用虛擬試驗代替真實試驗，并需要預先確定一套表達不同庫存期產品真實變化規律的實體模型和物理模型——本文稱之為真實模型.真實模型只用于產生虛擬試驗數據，文中不予展示.

虛擬試驗（下文簡稱試驗）的作法是，首先結合真實模型確定模型參數的設計值，然后在設計值上附加一個虛擬隨機量（以模仿真實的加工制造隨機性），將模型參數轉變為隨機變量，接著對該隨機變量進行抽樣，把抽樣結果視作試驗樣品真實狀態的模型參數，并據此求出各試驗樣品所對應的系統級輸出量的精確解.為了模仿真實試驗的測試隨機性，也通過附加虛擬隨機量的方法將精確解轉變為隨機性系統級輸出量，對該輸出量進行抽樣，便可獲取相當于真實試驗的測試結果.

表1 虛擬隨機變量的均方差Table 1 Deviations of virtual stochastic variables

4.3 模型校準

本例中需要校準的參數有兩類，其一是計算參數，其二是物理參數.參照驗證和確認的順序要求，先進行計算參數的校準，后進行物理參數的校準.因為新鮮材料和新鮮產品的試驗數據比較豐富，擬基于新鮮產品的試驗數據來校準計算參數，即首先根據工程任務對數值計算的準確性要求以及所耗費計算時間的承受能力，選擇恰當的空間步長，然后將時間步長和人為粘性系數確定為校準量，以新鮮產品的試驗結果為參照解，完成對時間步長和人為粘性系數的校準.物理參數的校準主要是對庫存老化后材料參數e1和γ2的校準，需要用到庫存產品的試驗數據.校準后的計算參數和物理參數形成一組固定搭配，用于新對象的數值模擬預測.

按以下步驟開展數值模擬校準：

1）對計算方法、計算參數、物理參數和可調參數進行應用需求分析，確定校準量以及參照解的獲取途徑.本例只在計算參數和物理參數中選取校準量（選取結果見下文），并假定分析對象是復雜系統，參照解來自試驗數據；

3）通過設計狀態下模擬結果和5個試驗結果的對比分析，得到優化的計算參數校準結果：初始網格寬度Δx＝0.1mm，時間步長Δt＝0.0016μs，人為粘性系數α＝1.5和β＝0.06，該粘性系數所對應的人為粘性模型為

其中q為粘性壓力，l為網格尺度，c為當地音速，Δ·v為速度的散度．計算參數校準以后，給出對應設計狀態下的模擬結果Dm＆s，結果列于表2；

表2 試驗結果和計算參數校準后的計算結果Table 2 Test data and M＆S result after com puting parameters calibrated

4）基于材料的老化模型以及庫存期為10年、30年和50年的系統級試驗與數值模擬的結果對比，校準物理模型參數.

老化模型給出的物理參數變化規律為

5）固化計算方法和各種校準參數的取值，完成數值模擬校準.

4.4 模型確認

在模型空間中，將庫存時間從0年～50年的區域定義為確認域，確認試驗分別對應0年、10年、20年、30年、40年和50年的庫存時間，每個庫存時間分別安排5次重復性試驗，但只對應一個數值模擬結果.

依據公式（7）進行確認環節的不確定度量化，當置信度β＝0.95時，t（1－β）/2，ν＝t0.025，4＝2.7764，確認域內的不確定度量化結果見表3.

表3 確認域內的數值模擬不確定度Table 3 M＆S uncertainties in validation domain

4.5 應用域內不確定度的量化

在模型空間中，將庫存時間大于50年的區域定義為應用域，該域內沒有新的系統級試驗.

為了量化式（3）的右端第一項，需要得到確認域內不確定度隨庫存時間的變化規律，本例用二次多項式來表達這個規律

根據表3中的信息，有

該超定方程的最小二范數解為

（a0，a1，a2）＝（1.415 7×10－1，2.6012×10－3，－1.3807×10－5）．

為量化式（3）的右端第二項，表4和表5分別列出了基于區間表達的不同庫存期模型參數的不確定度以及系統響應量D對模型參數的偏導數（即敏感系數），其中偶然不確定度源自加工隨機性導致的產品個體模型參數的隨機變化，而認知不確定度源自建模者對模型參數母體均值的認知缺陷.由于80年庫存期的材料級物理建模試驗較少，又沒有系統級試驗用于模型校準，因此在建模方面比0年～50年庫存期產品存在更多的認知不確定度，其中e1U增加了0.2MJ·kg－1，γ2U增加了0.2，即用式（8）－（10）求得的80年庫存期模型參數，對模型參數的真值存在更大的偏離.

表4 模型參數的不確定度Table 4 Uncertainties ofmodel parameters

表5 系統響應量對模型參數的偏導數Table 5 Partial derivatives of system response to model parameters

對庫存80年的產品，由模型參數不確定度增量傳遞到系統響應量上的不確定度為

根據式（3），對應80年庫存期的數值模擬不確定度為

假如將模型參數的認知不確定度整體傳遞到模擬結果中，則有

如果進一步將偶然不確定度也考慮在內，則僅由傳遞分析得到的數值模擬不確定度就高達0.92828mm，而且該結果還不能視為總的數值模擬不確定度.從這里可以看出，式（3）～式（7）所示的方法能有效利用試驗校準對系統級認知不確定度的消減作用，過濾來自試驗的偶然不確定度對量化結果的影響，合理地反映新對象數值模擬預測的不確定度.

4.6 不確定度量化方法的檢驗

為了檢驗上述方法的有效性，特意安排了應用域內對應80年庫存期的試驗.檢驗方法是，根據數值模擬和試驗數據，直接由式（7）所示的對比方法對不確定度進行評定，如果式（13）的量化結果大于該評定結果且大得不多，則在一定程度上說明量化結果滿足真值覆蓋和最小化原則，量化方法正確有效.

試驗方法同前，即根據真實模型確定80年庫存期的設計參數，然后結合精確解和抽樣方法給出5個重復性試驗結果.數值模擬繼續用式（10）所示的校準后的模型參數，并且數值模擬中模型參數的認知不確定度與虛擬實驗中模型參數的偶然不確定度符合表4，即80年庫存期與0～50年庫存期相比，數值模擬所用的e1和γ2對真實模型的偏離分別增加了0.2MJ·kg－1和0.2，這兩個增加的偏離在影響效果上均使D變大，能夠很好地反映不確定度的疊加效應，突出檢驗效果.

表6為庫存80年產品基于試驗對比的數值模擬不確定度評定結果，其中數值模擬所采用的模型參數為ρ1＝2 500kg·m－3，ρ2＝20 000kg·m－3，e1＝5.4MJ·kg－1，γ1＝3.0，γ2＝6.02，e1和γ2的取值來自式（8）－（10），不確定度的評定仍采用式（7）.

表6 基于試驗對比的不確定度評定Table 6 Uncertainty assessment based on com parison w ith test data

由表6知，試驗前用本文方法預測的數值模擬不確定度（0.381 65mm）稍大于試驗后由對比分析得到的不確定度評定結果（0.330 502mm），說明該例中的不確定度量化是成功的.

5 結論

當一個新對象不能進行實物試驗時，數值模擬預測及其不確定度量化結果就成為產品認證和評估的重要信息來源，而不確定度的傳遞分析往往被用作不確定度量化的常用手段.系統層級劃分越細，輸入的不確定因素也越多，僅靠傳遞分析則難以給出有實用價值的數值模擬不確定度.本文根據新對象的原型往往有一定試驗數據的工程現實，利用數值模擬的試驗校準能夠大量消減認知不確定度的作用機制，提出了在沒有新對象系統級試驗的條件下綜合利用校準后確認域對比信息和新增不確定度傳遞信息的數值模擬不確定度量化方法，案例表明該方法能夠較好地遵循不確定度量化的真值覆蓋原則和最小化原則，適用于新對象數值模擬預測的不確定度量化.

［1］ OberkampfW L，Roy C J.Verification and validation in science computing［M］.Cambridge University Press，2010：1－767.

［2］ Roy C J，OberkampfW L.A comprehensive framework for verification，validation，and uncertainty quantification in scientific computing［J］.Computer Methods in Applied Mechanics and Engineering，2011，200：2131－2144.

［3］ Ma Zhibo，Ying Yangjun，Zhu Jianshi.QMU certifyingmethod and its implementation［J］.Chinese Journalof Nuclear Science and Engineering，2009，29（1）：1－9.

［4］ Ma Zhibo，Li Haijie，Yin Jianwei，Huang Wenbin.Uncertainty quantification of numerical simulation for reliability analysis ［J］.Chinese Journal of Computational Physics，2014，31（4）：424－430.

［5］ Chen W，Baghdasaryan L，Buranathiti T，Cao J.Model validation via uncertainty propagation［J］.AIAA Journal，2004，42 （7）：1406－1415.

［6］ Helton J C.Conceptual and computational basis for the quantification of Margins and uncertaity［A］.Sandia Report，SAND2009－3055.

［7］ Liu Quan，Wang Ruili，Lin Zhong，Wen Wanzhi.Uncertainty quantification for JWL EOS parameters in explosive numerical simulation［J］.Chinese Journal of Explosion and Shock Waves，2013，33（6）：647－654.

［8］ Ma Zhibo，Zheng Miao，Yin Janwei，et al.Quantification of uncertainties in detonation simulations［J］.Chinese Journal of Computational Physics，2011，28（1）：66－74.

［9］ Higdon D，Kennedy M，et al.Combining field observation observations and simulation for calibration and prediction［J］. SIAM Journal of Scientific Computing，2004，26：448－466.

［10］ Helton J C.Uncertainty and sensitivity analysis in the presence of stochastic and subjective uncertainty［J］.Journal of Statistical Computation and Simulation，1997，57：3－76.

［11］ Babuska I，Nobile F，Tempone R.A systematic approach tomodel validation based on bayesian updates and prediction related rejection criteria［J］.Computer Methods in Applied Mechanics and Engineering，2008，197（29/32）：2517－2539.

Uncertainty Quantification of Numerical Simulations Subjected to Calibration

MA Zhibo1，YIN Jianwei1，LIHaijie2，LIU Quan1
（1.Institute of Applied Physics and Computational Mathematics，Beijing 100088，China；
2.Electromechanical Product Research and Development Center，Changsha 410100，China）

Propagation analysis often overrate uncertainty of numerical simulation，especially asmany uncertain inputs exist.Taking advantage of the reality that calibration can reduce epistemic uncertainty of system?level numerical simulation，amethod for uncertainty quantification is offered synthetically using comparison information based on available testdata and additional propagation information of modeling and simulation.An example with virtual test is displayed in which themethod is demonstrated and validated.

numerical simulation；uncertainty quantification；calibration；reliability certification

O242.1

A

2014－08－12；

2014－12－24

國家自然科學基金（11371066）、裝備預先研究（51319010207）及中國工程物理研究院科學技術發展基金（2012B0102010，2013A0101004）資助項目

馬智博（1963－），男，博士，研究員，主要從事計算物理和可靠性工程研究，E?mail：mazhibo＠iapcm.ac.cn

Received date： 2014－08－12；Revised date： 2014－12－24