校準(zhǔn)條件下的數(shù)值模擬不確定度量化方法

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

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100088；2.長(zhǎng)沙機(jī)電產(chǎn)品研究開(kāi)發(fā)中心，湖南長(zhǎng)沙 410100）

文章編號(hào)：1001?246X（2015）05?0514?09

校準(zhǔn)條件下的數(shù)值模擬不確定度量化方法

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

（1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100088；2.長(zhǎng)沙機(jī)電產(chǎn)品研究開(kāi)發(fā)中心，湖南長(zhǎng)沙 410100）

當(dāng)存在眾多不確定輸入因素時(shí)，不確定度的傳遞分析往往導(dǎo)致對(duì)數(shù)值模擬不確定度的過(guò)高估計(jì).利用校準(zhǔn)行為能夠消減系統(tǒng)級(jí)數(shù)值模擬中認(rèn)知不確定度的客觀機(jī)制，提出一個(gè)綜合利用已有系統(tǒng)級(jí)試驗(yàn)對(duì)比信息和新增建模與模擬傳遞信息的不確定度量化方法，結(jié)合一個(gè)虛擬試驗(yàn)的例子對(duì)該方法進(jìn)行展示和驗(yàn)證.

數(shù)值模擬；不確定度量化；校準(zhǔn)；可靠性認(rèn)證

0 引言

在投入應(yīng)用之前，建模與模擬（Modeling＆Simulation，M＆S）技術(shù)需要經(jīng)過(guò)驗(yàn)證（Verification）、確認(rèn)（Validation）和認(rèn)可（Accreditation）過(guò)程，以評(píng)估其可信度并獲取工程應(yīng)用所需的資質(zhì)［1－2］.事實(shí)上，即便經(jīng)過(guò)充分的驗(yàn)證和確認(rèn)，也難以完全發(fā)現(xiàn)并徹底解決M＆S中的所有問(wèn)題，再加上不可避免的離散誤差，對(duì)復(fù)雜物理過(guò)程的數(shù)值模擬往往會(huì)產(chǎn)生相對(duì)于真值的系統(tǒng)性偏離，因而常常通過(guò)校準(zhǔn)（Calibration）措施給予糾偏以提高模擬結(jié)果與試驗(yàn)結(jié)果的一致性.用經(jīng)過(guò)校準(zhǔn)的模擬程序?qū)π碌膹?fù)雜對(duì)象進(jìn)行模擬計(jì)算，是現(xiàn)實(shí)中普遍存在并將長(zhǎng)期存在的工程預(yù)測(cè)方法.當(dāng)不能進(jìn)行新對(duì)象的系統(tǒng)級(jí)試驗(yàn)時(shí)，對(duì)校準(zhǔn)后數(shù)值模擬的預(yù)測(cè)結(jié)果進(jìn)行不確定度量化，是可靠性認(rèn)證與評(píng)估所面臨的一個(gè)重要科學(xué)問(wèn)題［3］.

認(rèn)證的對(duì)象為抽象的事物類（如型號(hào)），評(píng)估的對(duì)象為具體的事物例（如產(chǎn)品）［4］.型號(hào)經(jīng)過(guò)認(rèn)證后獲得定型，定型產(chǎn)品經(jīng)過(guò)庫(kù)存和服役等過(guò)程后，可能發(fā)生狀態(tài)上的改變，往往需要分期評(píng)估.為全面考察各種工程因素的影響，新型號(hào)認(rèn)證一般會(huì)安排較多的系統(tǒng)級(jí)試驗(yàn)，試驗(yàn)數(shù)據(jù)可用于數(shù)值模擬的系統(tǒng)級(jí)校準(zhǔn)、確認(rèn)和不確定度量化.當(dāng)不能進(jìn)行新對(duì)象的系統(tǒng)級(jí)試驗(yàn)時(shí)，需要針對(duì)改變后的設(shè)計(jì)參數(shù)和新增加的工程因素進(jìn)行數(shù)值模擬及其模擬結(jié)果的不確定度量化.

如果不能進(jìn)行新對(duì)象的系統(tǒng)級(jí)試驗(yàn)，數(shù)值模擬不確定度量化的信息主要來(lái)自兩個(gè)方面：其一是對(duì)比信息，即通過(guò)原型對(duì)象的系統(tǒng)級(jí)試驗(yàn)數(shù)據(jù)與數(shù)值模擬的對(duì)比而得到的不確定度信息，利用模型空間中不確定度從確認(rèn)域到應(yīng)用域的外推技術(shù)，可以將該信息轉(zhuǎn)換為能夠直接反映新對(duì)象系統(tǒng)級(jí)數(shù)值模擬的不確定度信息［4］；其二是傳遞信息，即借助總體計(jì)算程序和擾動(dòng)分析方法，將較低層級(jí)建模與模擬不確定度信息傳遞到新對(duì)象系統(tǒng)級(jí)數(shù)值模擬的不確定度信息［5］.

可靠性分析往往受困于信息缺乏，多源信息的整合因此變得迫切.雖然對(duì)比信息和傳遞信息分別來(lái)自不同的認(rèn)知途徑，但反映的是同一個(gè)問(wèn)題，將二者進(jìn)行整合，既有現(xiàn)實(shí)需求，又有科學(xué)上的合理性和可行性.

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

數(shù)值模擬的不確定度主要為認(rèn)知不確定度，宜用區(qū)間理論來(lái)表達(dá)并進(jìn)行信息整合，不確定度量化應(yīng)當(dāng)遵循真值覆蓋原則和最小化原則［4］.如果量化結(jié)果僅利用外推后的對(duì)比信息，則會(huì)因?yàn)闆](méi)有充分考慮新對(duì)象建模和模擬所新增加的認(rèn)知缺陷，容易導(dǎo)致不確定度估計(jì)結(jié)果偏小而違背真值覆蓋原則，最終增加不合格對(duì)象被接收的風(fēng)險(xiǎn)，如果把外推后的對(duì)比信息和不確定度的整體傳遞信息直接相加或僅僅利用整體傳遞信息，均可能使不確定度估計(jì)結(jié)果被不合理地放大而違背最小化原則，從而增加合格對(duì)象被拒收的風(fēng)險(xiǎn).

根據(jù)認(rèn)知科學(xué)中對(duì)比信息和傳遞信息之間的容斥關(guān)系以及不確定度量化的基本原則，本文將經(jīng)過(guò)外推的系統(tǒng)級(jí)對(duì)比信息作為不確定度的基本組成部分，將新對(duì)象建模和模擬所新增加的不確定度傳遞到系統(tǒng)級(jí)的計(jì)算結(jié)果中，作為不確定度的增量部分，然后將二者進(jìn)行相加作為新對(duì)象數(shù)值模擬的不確定度，并結(jié)合一個(gè)黎曼問(wèn)題及其虛擬試驗(yàn)對(duì)該不確定度量化方法進(jìn)行展示和驗(yàn)證.

1 數(shù)值模擬校準(zhǔn)及其技術(shù)特征

提高計(jì)算結(jié)果與試驗(yàn)結(jié)果一致性的途徑有兩個(gè)，其一是全面提高認(rèn)知能力，降低數(shù)值模擬中每個(gè)技術(shù)環(huán)節(jié)的認(rèn)知不確定度，這也是數(shù)值模擬發(fā)展的理想途徑，其二是基于現(xiàn)有認(rèn)知能力，人為地讓不同技術(shù)環(huán)節(jié)所產(chǎn)生的認(rèn)知誤差相互補(bǔ)償.數(shù)值模擬校準(zhǔn)在較大程度上依賴于誤差補(bǔ)償機(jī)制，但這種補(bǔ)償作用只在較小的模型范圍內(nèi)有效，隨著模擬對(duì)象遠(yuǎn)離那些用以校準(zhǔn)的試驗(yàn)對(duì)象，數(shù)值模擬不確定度往往增加很快.

本文提及的模型包括實(shí)體模型和物理模型，前者指具體的特殊工程對(duì)象，包括部件的材料選用、密度與質(zhì)量、形狀與尺寸、動(dòng)作過(guò)程的初邊條件等工程設(shè)計(jì)參數(shù)，后者指抽象的一般物理規(guī)律，包括描述材料力學(xué)屬性的本構(gòu)關(guān)系和狀態(tài)方程、描述相關(guān)物理過(guò)程（如爆轟、湍流和傳熱等）的模型方程以及描述物質(zhì)運(yùn)動(dòng)普遍規(guī)律的三大守恒方程等.實(shí)體建模的不確定度往往決定于靜態(tài)測(cè)量的不確定度，物理建?？赡苌婕案嗟膭?dòng)態(tài)測(cè)量以及基于測(cè)量結(jié)果和近似理論的逆向運(yùn)算，因此，物理建模所引起的數(shù)值模擬不確定度一般較大.

結(jié)合試驗(yàn)進(jìn)行的校準(zhǔn)行為，主要體現(xiàn)在對(duì)模型形式、計(jì)算方法和選定的校準(zhǔn)參數(shù)進(jìn)行調(diào)整并把它們固定下來(lái)［2］，校準(zhǔn)參數(shù)一般來(lái)自計(jì)算參數(shù)、物理參數(shù)或可調(diào)參數(shù).可調(diào)參數(shù)可理解為缺乏明確數(shù)學(xué)和物理含義或缺乏確切賦值信息但對(duì)最終數(shù)值模擬結(jié)果有明顯影響的參數(shù).通過(guò)充分的驗(yàn)證和確認(rèn)，可有效減少可調(diào)參數(shù)的數(shù)量，提高數(shù)值模擬的預(yù)測(cè)能力，但由于數(shù)值計(jì)算對(duì)應(yīng)非零的時(shí)間和空間離散尺度，離散誤差不可避免，再加上模型形式的近似性和計(jì)算方法在模擬復(fù)雜問(wèn)題時(shí)能力的欠缺，可調(diào)參數(shù)的存在往往難以避免［9］.

校準(zhǔn)過(guò)程一般基于某一范圍內(nèi)的實(shí)體模型，本文將實(shí)體模型空間中的這個(gè)范圍稱為校準(zhǔn)域.校準(zhǔn)之后所進(jìn)行的確認(rèn)活動(dòng)，也有對(duì)應(yīng)的確認(rèn)域，其中模擬結(jié)果和試驗(yàn)結(jié)果的對(duì)比分析，被用于確認(rèn)域內(nèi)數(shù)值模擬的不確定度量化.

校準(zhǔn)過(guò)程完成以后，計(jì)算程序及其相關(guān)的模型參數(shù)即被固化，但固化通常具有相對(duì)性和階段性.相對(duì)性主要表現(xiàn)在描述物理規(guī)律的模型參數(shù)在取值上與特定計(jì)算方法和計(jì)算參數(shù)有關(guān)，并非完全基于對(duì)物理的認(rèn)識(shí)，而階段性表現(xiàn)在固化版本通常是暫時(shí)的，隨著數(shù)值模擬技術(shù)的發(fā)展，模型和參數(shù)也會(huì)隨之改變.

演繹和歸納是認(rèn)識(shí)事物的兩個(gè)基本方法，前者基于特定條件和普遍規(guī)律來(lái)推斷事物的發(fā)展結(jié)果，是不確定度傳遞分析的基礎(chǔ)，后者通過(guò)實(shí)踐和觀察來(lái)總結(jié)事物的發(fā)展結(jié)果，是不確定度對(duì)比分析的基礎(chǔ).對(duì)比信息直接來(lái)自系統(tǒng)層級(jí)的實(shí)踐活動(dòng)，因此可信度更高，在進(jìn)行信息整合時(shí)，對(duì)比信息的權(quán)重因子一般要遠(yuǎn)大于傳遞信息.

我們認(rèn)為，經(jīng)歷校準(zhǔn)后的數(shù)值模擬具有如下特點(diǎn)：

1）校準(zhǔn)域內(nèi)系統(tǒng)級(jí)數(shù)值模擬的不確定度能夠得到有效消減，但隨著實(shí)體模型遠(yuǎn)離校準(zhǔn)域，誤差補(bǔ)償作用逐漸減弱，新對(duì)象模擬結(jié)果對(duì)真值的背離會(huì)逐漸增大；

2）確認(rèn)域內(nèi)由對(duì)比分析得到的不確定度可外推至應(yīng)用域，但外推結(jié)果不能充分反映新對(duì)象建模和模擬所導(dǎo)致的新增不確定度；

3）系統(tǒng)級(jí)以下物理建模和數(shù)值計(jì)算的不確定度并不因?yàn)樾?zhǔn)而得到必然的消減，由傳遞分析得到的不確定度往往較大；

4）如果對(duì)比信息和傳遞信息均反映同一數(shù)值模擬的不確定度，進(jìn)行信息整合時(shí)，權(quán)重分配應(yīng)以對(duì)比信息為主；

5）當(dāng)確認(rèn)域外的新對(duì)象不能進(jìn)行試驗(yàn)時(shí)，數(shù)值模擬不確定度存在兩個(gè)獨(dú)立的信息來(lái)源：其一是系統(tǒng)級(jí)確認(rèn)域內(nèi)來(lái)自對(duì)比的不確定度信息，其二是系統(tǒng)級(jí)以下新增的建模和模擬不確定度所傳遞到系統(tǒng)級(jí)的不確定度信息，二者具有可加性，可基于區(qū)間理論進(jìn)行信息整合.

2 校準(zhǔn)條件下的不確定度量化

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

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

當(dāng)沒(méi)有新對(duì)象的系統(tǒng)級(jí)試驗(yàn)時(shí)，無(wú)直接的對(duì)比信息可利用，不確定度量化的可能方案有：

1）僅僅利用傳遞信息；

2）僅僅利用外推后的對(duì)比信息；

3）同時(shí)利用以上兩種信息.

第一方案未利用原有的系統(tǒng)級(jí)對(duì)比信息，由傳遞給出的認(rèn)知不確定度不能得到有效消減，傳遞過(guò)程產(chǎn)生的不確定度難以量化，低層級(jí)數(shù)值計(jì)算造成的不確定度也難以傳遞到系統(tǒng)層級(jí).第二方案可能漏掉新對(duì)象建模和模擬所新產(chǎn)生的不確定度.第三方案有更加合理的思想基礎(chǔ)，但需要設(shè)計(jì)一套恰當(dāng)?shù)膶?shí)施方法，最大限度地減少上述問(wèn)題造成的影響.

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

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

其中1Upropagation，2Upropagation，3Upropagation分別表示實(shí)體建模、物理建模和數(shù)值計(jì)算三個(gè)技術(shù)環(huán)節(jié)所傳遞到系統(tǒng)級(jí)的數(shù)值模擬不確定度.

新對(duì)象數(shù)

綜上所述，可據(jù)如下步驟進(jìn)行數(shù)值模擬不確定度的量化：

1）結(jié)合校準(zhǔn)試驗(yàn)，進(jìn)行數(shù)值模擬校準(zhǔn)，然后完成對(duì)計(jì)算程序和相關(guān)參數(shù)的固化；

2）結(jié)合確認(rèn)試驗(yàn)，進(jìn)行系統(tǒng)級(jí)確認(rèn)域內(nèi)的對(duì)比分析和不確定度量化；

3）根據(jù)不確定度隨實(shí)體模型參數(shù)的變化規(guī)律，用確認(rèn)域內(nèi)的不確定度信息外推新對(duì)象數(shù)值模擬的不確定度；

5）根據(jù)公式（3）得到新對(duì)象數(shù)值模擬的總體不確定度.

3 基于對(duì)比分析的不確定度量化

試驗(yàn)對(duì)象的加工隨機(jī)性和測(cè)試隨機(jī)性都會(huì)使試驗(yàn)結(jié)果產(chǎn)生偶然不確定度，但該不確定度并非由數(shù)值模擬產(chǎn)生，因此，基于對(duì)比分析的不確定度量化方法，要恰當(dāng)?shù)貥?gòu)造統(tǒng)計(jì)量，屏除這類不確定度的影響，同時(shí)有效地揭示真正的數(shù)值模擬不確定度信息.

在工程實(shí)踐中，基于試驗(yàn)的對(duì)比分析有兩種情形：

1）一對(duì)一，即一次數(shù)值模擬對(duì)應(yīng)一次試驗(yàn)；

2）一對(duì)多，即一次數(shù)值模擬對(duì)應(yīng)多次試驗(yàn).

在一對(duì)一情形中，需要對(duì)每次試驗(yàn)的對(duì)象進(jìn)行單獨(dú)測(cè)量，然后將測(cè)量結(jié)果用于實(shí)體建模和物理建模，每個(gè)試驗(yàn)對(duì)象均有相應(yīng)的數(shù)值模擬結(jié)果和試驗(yàn)結(jié)果，在扣除測(cè)試誤差后，兩者之間的差別能夠直接地反映試驗(yàn)對(duì)象的數(shù)值模擬誤差.如果試驗(yàn)和模擬次數(shù)均為n，可按下式提取數(shù)值模擬的不確定度：

一對(duì)多情形往往對(duì)應(yīng)于重復(fù)性試驗(yàn)，僅基于設(shè)計(jì)方案的公稱參數(shù)或?qū)嶋H加工的中值參數(shù)進(jìn)行實(shí)體建模和物理建模，最終只產(chǎn)生一套數(shù)值模擬結(jié)果.雖然每一個(gè)試驗(yàn)對(duì)象都按照同樣的設(shè)計(jì)方案進(jìn)行加工制造，但產(chǎn)品加工和試驗(yàn)測(cè)試的隨機(jī)性會(huì)導(dǎo)致試驗(yàn)結(jié)果的隨機(jī)性.為了屏除隨機(jī)性干擾信息，本文建議從數(shù)值模擬結(jié)果與試驗(yàn)均值之間的差異中挖掘數(shù)值模擬的不確定度信息，具體公式

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

當(dāng)所有試驗(yàn)對(duì)象均處于某設(shè)計(jì)點(diǎn)附近且n足夠大時(shí)，式（4）、（7）給出的結(jié)果趨于一致.

4 實(shí)例演示

4.1 問(wèn)題描述

圖1所示的實(shí)體模型對(duì)應(yīng)一維黎曼問(wèn)題，共有兩個(gè)物質(zhì)區(qū)域，左端區(qū)域相當(dāng)于高壓氣體，其初始時(shí)刻的壓強(qiáng)、密度、內(nèi)能和速度分別為p1，ρ1，e1，v1，右端區(qū)域相當(dāng)于重金屬材料，其初始時(shí)刻的壓強(qiáng)、密度、內(nèi)能和速度分別為p2，ρ2，e2，v2．兩區(qū)在零時(shí)刻的初速度均為零，從零時(shí)刻開(kāi)始，在物質(zhì)界面處產(chǎn)生一個(gè)向左傳播的膨脹波和向右傳播的沖擊波，金屬材料在高壓氣體的推動(dòng)下向右運(yùn)動(dòng)，本例所關(guān)心的系統(tǒng)級(jí)輸出量是10μs時(shí)物質(zhì)界面移動(dòng)的距離D.

圖1 實(shí)體模型示意圖Fig.1 Sketchmap of entitymodel

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

本例將上述實(shí)體模型看作一個(gè)產(chǎn)品型號(hào)，并假定庫(kù)存老化僅導(dǎo)致兩個(gè)物理參數(shù)發(fā)生改變，其一是氣體的初始內(nèi)能e1（初始?jí)簭?qiáng)與初始內(nèi)能相關(guān)，不作為獨(dú)立的模型變量），其二是金屬的多方指數(shù)γ2．因二者均為庫(kù)存時(shí)間t的函數(shù)，可基于庫(kù)存時(shí)間建立實(shí)體模型的模型空間，并在該模型空間中定義校準(zhǔn)域、確認(rèn)域和應(yīng)用域.

要解決的問(wèn)題是，當(dāng)不能對(duì)庫(kù)存產(chǎn)品進(jìn)行新的系統(tǒng)級(jí)試驗(yàn)時(shí)，如何基于已有的0年～50年庫(kù)存期產(chǎn)品的系統(tǒng)級(jí)試驗(yàn)結(jié)果及其對(duì)應(yīng)的數(shù)值模擬結(jié)果，量化80年庫(kù)存期產(chǎn)品數(shù)值模擬預(yù)測(cè)的不確定度.

4.2 試驗(yàn)設(shè)計(jì)

模型校準(zhǔn)、模型確認(rèn)以及對(duì)不確定度量化方法的檢驗(yàn)將涉及試驗(yàn)，該例存在系統(tǒng)級(jí)輸出量的精確解D?，因此擬用虛擬試驗(yàn)代替真實(shí)試驗(yàn)，并需要預(yù)先確定一套表達(dá)不同庫(kù)存期產(chǎn)品真實(shí)變化規(guī)律的實(shí)體模型和物理模型——本文稱之為真實(shí)模型.真實(shí)模型只用于產(chǎn)生虛擬試驗(yàn)數(shù)據(jù)，文中不予展示.

虛擬試驗(yàn)（下文簡(jiǎn)稱試驗(yàn)）的作法是，首先結(jié)合真實(shí)模型確定模型參數(shù)的設(shè)計(jì)值，然后在設(shè)計(jì)值上附加一個(gè)虛擬隨機(jī)量（以模仿真實(shí)的加工制造隨機(jī)性），將模型參數(shù)轉(zhuǎn)變?yōu)殡S機(jī)變量，接著對(duì)該隨機(jī)變量進(jìn)行抽樣，把抽樣結(jié)果視作試驗(yàn)樣品真實(shí)狀態(tài)的模型參數(shù)，并據(jù)此求出各試驗(yàn)樣品所對(duì)應(yīng)的系統(tǒng)級(jí)輸出量的精確解.為了模仿真實(shí)試驗(yàn)的測(cè)試隨機(jī)性，也通過(guò)附加虛擬隨機(jī)量的方法將精確解轉(zhuǎn)變?yōu)殡S機(jī)性系統(tǒng)級(jí)輸出量，對(duì)該輸出量進(jìn)行抽樣，便可獲取相當(dāng)于真實(shí)試驗(yàn)的測(cè)試結(jié)果.

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

4.3 模型校準(zhǔn)

本例中需要校準(zhǔn)的參數(shù)有兩類，其一是計(jì)算參數(shù)，其二是物理參數(shù).參照驗(yàn)證和確認(rèn)的順序要求，先進(jìn)行計(jì)算參數(shù)的校準(zhǔn)，后進(jìn)行物理參數(shù)的校準(zhǔn).因?yàn)樾迈r材料和新鮮產(chǎn)品的試驗(yàn)數(shù)據(jù)比較豐富，擬基于新鮮產(chǎn)品的試驗(yàn)數(shù)據(jù)來(lái)校準(zhǔn)計(jì)算參數(shù)，即首先根據(jù)工程任務(wù)對(duì)數(shù)值計(jì)算的準(zhǔn)確性要求以及所耗費(fèi)計(jì)算時(shí)間的承受能力，選擇恰當(dāng)?shù)目臻g步長(zhǎng)，然后將時(shí)間步長(zhǎng)和人為粘性系數(shù)確定為校準(zhǔn)量，以新鮮產(chǎn)品的試驗(yàn)結(jié)果為參照解，完成對(duì)時(shí)間步長(zhǎng)和人為粘性系數(shù)的校準(zhǔn).物理參數(shù)的校準(zhǔn)主要是對(duì)庫(kù)存老化后材料參數(shù)e1和γ2的校準(zhǔn)，需要用到庫(kù)存產(chǎn)品的試驗(yàn)數(shù)據(jù).校準(zhǔn)后的計(jì)算參數(shù)和物理參數(shù)形成一組固定搭配，用于新對(duì)象的數(shù)值模擬預(yù)測(cè).

按以下步驟開(kāi)展數(shù)值模擬校準(zhǔn)：

1）對(duì)計(jì)算方法、計(jì)算參數(shù)、物理參數(shù)和可調(diào)參數(shù)進(jìn)行應(yīng)用需求分析，確定校準(zhǔn)量以及參照解的獲取途徑.本例只在計(jì)算參數(shù)和物理參數(shù)中選取校準(zhǔn)量（選取結(jié)果見(jiàn)下文），并假定分析對(duì)象是復(fù)雜系統(tǒng)，參照解來(lái)自試驗(yàn)數(shù)據(jù)；

3）通過(guò)設(shè)計(jì)狀態(tài)下模擬結(jié)果和5個(gè)試驗(yàn)結(jié)果的對(duì)比分析，得到優(yōu)化的計(jì)算參數(shù)校準(zhǔn)結(jié)果：初始網(wǎng)格寬度Δx＝0.1mm，時(shí)間步長(zhǎng)Δt＝0.0016μs，人為粘性系數(shù)α＝1.5和β＝0.06，該粘性系數(shù)所對(duì)應(yīng)的人為粘性模型為

其中q為粘性壓力，l為網(wǎng)格尺度，c為當(dāng)?shù)匾羲?，Δ·v為速度的散度．計(jì)算參數(shù)校準(zhǔn)以后，給出對(duì)應(yīng)設(shè)計(jì)狀態(tài)下的模擬結(jié)果Dm＆s，結(jié)果列于表2；

表2 試驗(yàn)結(jié)果和計(jì)算參數(shù)校準(zhǔn)后的計(jì)算結(jié)果Table 2 Test data and M＆S result after com puting parameters calibrated

4）基于材料的老化模型以及庫(kù)存期為10年、30年和50年的系統(tǒng)級(jí)試驗(yàn)與數(shù)值模擬的結(jié)果對(duì)比，校準(zhǔn)物理模型參數(shù).

老化模型給出的物理參數(shù)變化規(guī)律為

5）固化計(jì)算方法和各種校準(zhǔn)參數(shù)的取值，完成數(shù)值模擬校準(zhǔn).

4.4 模型確認(rèn)

在模型空間中，將庫(kù)存時(shí)間從0年～50年的區(qū)域定義為確認(rèn)域，確認(rèn)試驗(yàn)分別對(duì)應(yīng)0年、10年、20年、30年、40年和50年的庫(kù)存時(shí)間，每個(gè)庫(kù)存時(shí)間分別安排5次重復(fù)性試驗(yàn)，但只對(duì)應(yīng)一個(gè)數(shù)值模擬結(jié)果.

依據(jù)公式（7）進(jìn)行確認(rèn)環(huán)節(jié)的不確定度量化，當(dāng)置信度β＝0.95時(shí)，t（1－β）/2，ν＝t0.025，4＝2.7764，確認(rèn)域內(nèi)的不確定度量化結(jié)果見(jiàn)表3.

表3 確認(rèn)域內(nèi)的數(shù)值模擬不確定度Table 3 M＆S uncertainties in validation domain

4.5 應(yīng)用域內(nèi)不確定度的量化

在模型空間中，將庫(kù)存時(shí)間大于50年的區(qū)域定義為應(yīng)用域，該域內(nèi)沒(méi)有新的系統(tǒng)級(jí)試驗(yàn).

為了量化式（3）的右端第一項(xiàng)，需要得到確認(rèn)域內(nèi)不確定度隨庫(kù)存時(shí)間的變化規(guī)律，本例用二次多項(xiàng)式來(lái)表達(dá)這個(gè)規(guī)律

根據(jù)表3中的信息，有

該超定方程的最小二范數(shù)解為

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

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

表4 模型參數(shù)的不確定度Table 4 Uncertainties ofmodel parameters

表5 系統(tǒng)響應(yīng)量對(duì)模型參數(shù)的偏導(dǎo)數(shù)Table 5 Partial derivatives of system response to model parameters

對(duì)庫(kù)存80年的產(chǎn)品，由模型參數(shù)不確定度增量傳遞到系統(tǒng)響應(yīng)量上的不確定度為

根據(jù)式（3），對(duì)應(yīng)80年庫(kù)存期的數(shù)值模擬不確定度為

假如將模型參數(shù)的認(rèn)知不確定度整體傳遞到模擬結(jié)果中，則有

如果進(jìn)一步將偶然不確定度也考慮在內(nèi)，則僅由傳遞分析得到的數(shù)值模擬不確定度就高達(dá)0.92828mm，而且該結(jié)果還不能視為總的數(shù)值模擬不確定度.從這里可以看出，式（3）～式（7）所示的方法能有效利用試驗(yàn)校準(zhǔn)對(duì)系統(tǒng)級(jí)認(rèn)知不確定度的消減作用，過(guò)濾來(lái)自試驗(yàn)的偶然不確定度對(duì)量化結(jié)果的影響，合理地反映新對(duì)象數(shù)值模擬預(yù)測(cè)的不確定度.

4.6 不確定度量化方法的檢驗(yàn)

為了檢驗(yàn)上述方法的有效性，特意安排了應(yīng)用域內(nèi)對(duì)應(yīng)80年庫(kù)存期的試驗(yàn).檢驗(yàn)方法是，根據(jù)數(shù)值模擬和試驗(yàn)數(shù)據(jù)，直接由式（7）所示的對(duì)比方法對(duì)不確定度進(jìn)行評(píng)定，如果式（13）的量化結(jié)果大于該評(píng)定結(jié)果且大得不多，則在一定程度上說(shuō)明量化結(jié)果滿足真值覆蓋和最小化原則，量化方法正確有效.

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

表6為庫(kù)存80年產(chǎn)品基于試驗(yàn)對(duì)比的數(shù)值模擬不確定度評(píng)定結(jié)果，其中數(shù)值模擬所采用的模型參數(shù)為ρ1＝2 500kg·m－3，ρ2＝20 000kg·m－3，e1＝5.4MJ·kg－1，γ1＝3.0，γ2＝6.02，e1和γ2的取值來(lái)自式（8）－（10），不確定度的評(píng)定仍采用式（7）.

表6 基于試驗(yàn)對(duì)比的不確定度評(píng)定Table 6 Uncertainty assessment based on com parison w ith test data

由表6知，試驗(yàn)前用本文方法預(yù)測(cè)的數(shù)值模擬不確定度（0.381 65mm）稍大于試驗(yàn)后由對(duì)比分析得到的不確定度評(píng)定結(jié)果（0.330 502mm），說(shuō)明該例中的不確定度量化是成功的.

5 結(jié)論

當(dāng)一個(gè)新對(duì)象不能進(jìn)行實(shí)物試驗(yàn)時(shí)，數(shù)值模擬預(yù)測(cè)及其不確定度量化結(jié)果就成為產(chǎn)品認(rèn)證和評(píng)估的重要信息來(lái)源，而不確定度的傳遞分析往往被用作不確定度量化的常用手段.系統(tǒng)層級(jí)劃分越細(xì)，輸入的不確定因素也越多，僅靠傳遞分析則難以給出有實(shí)用價(jià)值的數(shù)值模擬不確定度.本文根據(jù)新對(duì)象的原型往往有一定試驗(yàn)數(shù)據(jù)的工程現(xiàn)實(shí)，利用數(shù)值模擬的試驗(yàn)校準(zhǔn)能夠大量消減認(rèn)知不確定度的作用機(jī)制，提出了在沒(méi)有新對(duì)象系統(tǒng)級(jí)試驗(yàn)的條件下綜合利用校準(zhǔn)后確認(rèn)域?qū)Ρ刃畔⒑托略霾淮_定度傳遞信息的數(shù)值模擬不確定度量化方法，案例表明該方法能夠較好地遵循不確定度量化的真值覆蓋原則和最小化原則，適用于新對(duì)象數(shù)值模擬預(yù)測(cè)的不確定度量化.

［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

國(guó)家自然科學(xué)基金（11371066）、裝備預(yù)先研究（51319010207）及中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金（2012B0102010，2013A0101004）資助項(xiàng)目

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

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