梁 霄,王瑞利

(1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100094；2.山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，山東青島 266590)

1 引 言

基于建模與模擬(Modeling and Simulation,M&S)的設(shè)計(jì)可以明顯縮短產(chǎn)品研發(fā)周期、降低其生產(chǎn)成本，在生產(chǎn)和研發(fā)過程中起著重要的作用，尤其是試驗(yàn)成本很高、風(fēng)險(xiǎn)很大和危害人身安全的研究領(lǐng)域。如何提高建模的逼真度及模擬結(jié)果的精度，從而通過物理或者數(shù)學(xué)模型模擬客觀世界及其演變過程，已成為一個(gè)亟待解決的重要問題[1-8]。

研究者總是希望M&S方法能夠?qū)ρ芯靠腕w進(jìn)行逼真的刻畫和精確的描述，然而，實(shí)際研究過程中，誤差和不確定度的存在是不可避免的[8]。 特別是對(duì)于一些復(fù)雜的非線性系統(tǒng)，初值或邊值的微小差異有時(shí)會(huì)引起模擬結(jié)果劇烈的變化，影響模型的預(yù)測(cè)效果。M&S方法的關(guān)鍵任務(wù)之一就在于掌握誤差與不確定度的規(guī)律，通過采用比較完善的M&S技術(shù)，使用更加可靠的算法，將不確定度和誤差減小到最低限度。為了達(dá)到這一目的，分析M&S不確定度的來源及其性質(zhì)是非常重要的。了解不確定度對(duì)計(jì)算結(jié)果的影響，可以確定最終結(jié)果的可信度。

不確定度量化(UQ)是一門蓬勃發(fā)展的學(xué)科。美國(guó)非常重視UQ技術(shù)的研究，近幾年一直將UQ技術(shù)作為數(shù)值模擬軟件及其應(yīng)用研究的重點(diǎn)。1991年暫停地下核試驗(yàn)后，美國(guó)為了維持武器庫的水平與能力，保證武器的可靠性，開展了先進(jìn)模擬和計(jì)算計(jì)劃，進(jìn)一步刺激了UQ的發(fā)展。另外，2001年隸屬于美國(guó)能源部的三大國(guó)家核武器安全實(shí)驗(yàn)室(SNL,LANL,LLNL)共同倡導(dǎo)將UQ方法應(yīng)用于核武器庫存安全管理，給出了認(rèn)知不確定度和偶然不確定度在反應(yīng)堆故障排除和放射性核廢料處理中的應(yīng)用案例[9-11]。本研究以爆炸波問題為例，利用非嵌入多項(xiàng)式混沌方法進(jìn)行偶然不確定度的量化，并將結(jié)果與Monte Carlo(MC)所得結(jié)果比較，以期為非嵌入多項(xiàng)式系統(tǒng)不確定度的量化提供參考。

2 不確定度的分類和量化方法

如圖1所示，不確定度一般分為偶然不確定度和認(rèn)知不確定度。偶然不確定度是由于事物固有的隨機(jī)性而引起的，有時(shí)也稱為可變性和隨機(jī)不確定性，是系統(tǒng)固有的特點(diǎn)，即使在理論上也不能約減，如火箭助推器系統(tǒng)中質(zhì)量特性和幾何形狀的變異。目前對(duì)偶然不確定度的數(shù)學(xué)描述是采用概率論的方法，給出輸出量的概率描述。增加數(shù)據(jù)的樣本容量有助于更準(zhǔn)確地確定樣本的分布類型，但并不會(huì)減小不確定度。認(rèn)知不確定度是由于缺乏知識(shí)導(dǎo)致的不確定度，有時(shí)也稱為知識(shí)的不確定性、模型形式的不確定性或主觀不確定性。認(rèn)知不確定度分兩類：一類是認(rèn)可不確定度，即所謂“知道的不知道”，來源于對(duì)物理現(xiàn)象的理解較差，如斷裂力學(xué)； 另一類是盲目不確定度，即所謂“不知道的不知道”，例如軟件中的程序代碼出錯(cuò)[8]。

圖1 不確定度分類Fig.1 Classification of uncertainty

不確定度量化方法是利用各種數(shù)學(xué)和統(tǒng)計(jì)學(xué)方法，對(duì)隨機(jī)因素或認(rèn)知缺陷盡可能地給出分布或區(qū)間描述，建立M&S中偶然或認(rèn)知不確定性因素與系統(tǒng)響應(yīng)量之間的關(guān)系，進(jìn)而量化影響系統(tǒng)性能的不確定性，給出系統(tǒng)響應(yīng)量性能分析的不確定度度量方法。

在計(jì)算流體力學(xué)和計(jì)算固體力學(xué)領(lǐng)域中，不確定度量化的概念和方法發(fā)展迅速[8]。許多方法都可以應(yīng)用于偶然不確定度的量化，其中敏感度分析方法、抽樣方法、代理模型方法、隨機(jī)微分方程方法、參數(shù)率定與優(yōu)化方法已經(jīng)應(yīng)用于復(fù)雜工程中。MC方法是出現(xiàn)最早、應(yīng)用最廣泛的一種抽樣方法，并且使用直接，不需要任何假設(shè)，可以得出模型在不確定意義下的精確解，通常作為基準(zhǔn)解[12]。但是MC方法的計(jì)算成本過高，還會(huì)遇到維數(shù)災(zāi)難問題，因此人們已經(jīng)開始尋找MC的替代方法，如隨機(jī)譜方法。

隨機(jī)譜方法是使用譜逼近隨機(jī)微分方程或參量，然后將其分解為獨(dú)立的確定性分量和隨機(jī)分量來量化不確定性。其數(shù)學(xué)思想是數(shù)學(xué)模型中的每一個(gè)參量利用正交多項(xiàng)式(如Hermite多項(xiàng)式)展開成無窮級(jí)數(shù)項(xiàng)，在實(shí)際應(yīng)用中取有限項(xiàng)[13-15]。根據(jù)求解系數(shù)方法的不同，將多項(xiàng)式混沌分為嵌入多項(xiàng)式混沌(Intrusive Polynomial Chaos,IPC)和非嵌入多項(xiàng)式混沌(Non-Intrusive Polynomial Chaos,NIPC)。IPC無法利用已有程序，需要對(duì)已有的數(shù)值求解程序進(jìn)行大量修改或者重新編制計(jì)算程序。

3 非嵌入多項(xiàng)式混沌方法

NIPC是把已有的數(shù)值求解程序作為一個(gè)黑匣子，在隨機(jī)空間里通過一定的抽樣方法獲得若干個(gè)樣本點(diǎn)，將樣本點(diǎn)輸入到確定性程序求解，然后對(duì)確定性輸出結(jié)果進(jìn)行統(tǒng)計(jì)分析，以獲得相關(guān)數(shù)值求解結(jié)果的統(tǒng)計(jì)特征，進(jìn)而評(píng)估輸入?yún)?shù)或計(jì)算條件的不確定性在計(jì)算過程中傳播的影響。NIPC采用已有的數(shù)值模擬程序，不需要對(duì)控制方程進(jìn)行修改，不需要重新編寫程序。

NIPC不確定度評(píng)估的計(jì)算流程如下。

(1) 選取參數(shù)λ服從何種分布。

(2) 利用相應(yīng)標(biāo)準(zhǔn)隨機(jī)變量，確定參數(shù)λ多項(xiàng)式混沌譜展開，如λ(ξ)=λ0+sξ。

(3) 對(duì)參數(shù)構(gòu)成的隨機(jī)變量ξ進(jìn)行抽樣{θ1，θ2，…，θM+1}，將其帶入?yún)?shù)λ多項(xiàng)式混沌譜展開式中，得到參數(shù)λ的抽樣{ω1，ω2，…，ωM+1}，對(duì)每一次抽樣，模型參數(shù)值為確定值ωj，j=1,…，M+1。

(4) 對(duì)每個(gè)抽樣確定后的模型參數(shù)值，利用確定性程序(例如CFD程序)計(jì)算得到相應(yīng)的數(shù)值解u(x,t)，如第j次抽樣得到解為[u(x,t)]j。

(1)式可采用眾多積分公式求出，例如：假定正交多項(xiàng)式歸一，采用Latin等概率抽樣，則

(2)

4 偶然不確定度量化在爆炸波問題中的應(yīng)用

4.1 爆炸波問題的控制方程和迎風(fēng)格式

一維Lagrange坐標(biāo)系中爆炸波問題的控制方程為[16]

(3)

(3)式也可以寫成非守恒形式

(4)

4.2 爆炸波初邊值問題的不確定度量化

與Lax和Sod問題類似，爆炸波問題也是一種激波管問題。根據(jù)炸藥PBX9401的特性，將初始狀態(tài)等效于炸藥爆轟波CJ狀態(tài)，建立一維平面模型：在無限長(zhǎng)管道中，以x=0為界，左側(cè)是炸藥，采用高壓氣體模型描述，右側(cè)是低壓氣體。初始時(shí)刻，兩側(cè)物質(zhì)均處于靜止?fàn)顟B(tài)，中間用薄膜隔離。t=0時(shí)刻，薄膜發(fā)生破裂，產(chǎn)生左行稀疏波和右行激波。假設(shè)該爆炸波問題的計(jì)算區(qū)域?yàn)閇-1,1]，初值條件如下：在x≤0區(qū)域，炸藥密度ρL服從期望μ為1.4 g/cm3、標(biāo)準(zhǔn)差σ為5.0×10-4g/cm3的正態(tài)分布，速度uL=-2.6 km/s，壓強(qiáng)pL=37.18 GPa，多方指數(shù)γL=3.1；在x>0區(qū)域，氣體密度ρR=15.3 g/cm3，速度uR=0，壓強(qiáng)pR=1.0 GPa，多方指數(shù)γR=1.4。

該問題的特點(diǎn)是左、右兩側(cè)氣體的多方指數(shù)不同，屬于多介質(zhì)激波管問題，并且有解析解[17]。該問題的解析解表明，強(qiáng)間斷橫跨激波面和激波后有一個(gè)窄的密度峰，可以用于驗(yàn)證流體力學(xué)程序的計(jì)算能力和格式精度。對(duì)于左端炸藥密度，由于加工過程中炸藥顆粒凝結(jié)的不均勻性，其密度存在隨機(jī)性。假設(shè)炸藥密度的概率密度函數(shù)(PDF)為

(5)

圖2 炸藥密度分布直方圖Fig.2 Histogram of explosive density distribution

確定型爆炸波方程通常使用數(shù)學(xué)模型的單一解描述。對(duì)于不確定系統(tǒng)，這是不現(xiàn)實(shí)的，實(shí)際上需要執(zhí)行一系列運(yùn)算評(píng)估模型中輸入不確定度對(duì)輸出不確定度的影響，即不確定度傳播。確定型方程利用迎風(fēng)格式，計(jì)算到T=0.10時(shí)刻，得到1 000組狀態(tài)響應(yīng)量構(gòu)成的“馬尾圖”。計(jì)算結(jié)果的信息量巨大，必須對(duì)其量化以提取有用的信息。

圖3 1 000組抽樣下爆炸波問題的計(jì)算結(jié)果Fig.3 Results of blast wave problem solved with 1 000 times sampling by NIPC method

4.3 爆炸波問題的不確定度傳播

采用NIPC和MC方法對(duì)爆炸波初值問題開展不確定度的量化，得到T=0.10時(shí)刻計(jì)算區(qū)域中密度、壓力、速度、能量的期望和標(biāo)準(zhǔn)差如圖4所示，其中E為期望，σ為標(biāo)準(zhǔn)差。由圖4可知，MC方法和NIPC方法吻合較好，說明NIPC應(yīng)用于不確定度量化分析是可行的。

圖4 T=0.10時(shí)刻爆炸波問題中密度、壓強(qiáng)、速度和能量的期望與標(biāo)準(zhǔn)差Fig.4 Expectations and standard variances of density,pressure,velocity and internal energy in blast wave problem at T=0.10

4.4 激波位置的概率密度函數(shù)

為進(jìn)一步驗(yàn)證NIPC方法的有效性，利用NIPC方法研究初邊值輸入不確定度對(duì)爆炸波問題激波位置的影響，計(jì)算得到T=0.10時(shí)刻激波位置的概率密度函數(shù)，如圖5所示?？梢钥闯?，激波位置近似服從正態(tài)分布，激波位置的期望為0.19，標(biāo)準(zhǔn)差為4.13×10-4。圖6為采用MC方法取樣106次與采用NIPC取樣103次計(jì)算所得激波位置的概率密度函數(shù)，通過對(duì)比可以看出，MC方法的概率密度函數(shù)和NIPC方法得到的概率密度函數(shù)基本吻合，說明了NIPC方法的可靠性和有效性。

圖5 激波位置的概率密度分布函數(shù)Fig.5 Probability density funtion of shock wave position by NIPC at T=0.10

圖6 激波位置的概率密度分布函數(shù)Fig.6 Comparison of probability density funtion of shock wave position by MC and NIPC at T=0.10

5 結(jié) 論

給出建模中輸入初邊值偶然不確定度抽樣方法及量化結(jié)果，并利用MC和NIPC方法計(jì)算得到偶然不確定度對(duì)輸出結(jié)果期望、標(biāo)準(zhǔn)差及激波位置的影響，驗(yàn)證了NIPC方法在偶然不確定度量化方面的有效性。 此方法對(duì)復(fù)雜工程M&S中偶然不確定度量化具有重要的參考價(jià)值。

[1] 楊振虎.假設(shè)檢驗(yàn)與CFD不確定度量化分析 [J].航空計(jì)算技術(shù),2003,33(2):5-8.

YANG Z H.Hypothesis tests and analysis of CFD uncertainty quantification [J].Aeronautical Computer Technique,2003,33(2):5-8.

[2] 王瑞利,江 松.多物理耦合非線性偏微分方程與數(shù)值解不確定度量化數(shù)學(xué)方法 [J].中國(guó)科學(xué)A輯,2015,45(5):723-728.

WANG R L,JIANG S.Mathematical methods for uncertainty quantification of the nonlinear partial differential equation and numerical solution [J].Science China Series A,2015,45(5):723-728.

[3] 劉 全,王瑞利,林 忠,等.流體力學(xué)拉氏程序收斂性及數(shù)值計(jì)算不確定度初探 [J].計(jì)算物理,2013,30(3):346-352.

LIU Q,WANG R L,LIN Z,et al.Asymptotic convergence analysis and quantification of uncertainty in Lagrangian computation [J].Chinese Journal of Computational Physics,2013,30(3):346-352.

[4] 王瑞利,劉 全,溫萬治.非嵌入式多項(xiàng)式混沌法在爆轟產(chǎn)物JWL參數(shù)評(píng)估中的應(yīng)用 [J].爆炸與沖擊,2015,35(1):9-15.

WANG R L,LIU Q,WEN W Z.Uncertainty quantification for JWL EOS parameters using non-intrusive polynomial chaos [J].Explosion and Shock Waves,2015,35(1):9-15.

[5] 唐 煒,王玉明.復(fù)雜系統(tǒng)關(guān)鍵特征參數(shù)確定方法 [J].信息與電子工程,2011,9(1):83-86.

TANG W,WANG Y M.A method of confirming critical characteristic parameters for complex system [J].Information and Electrical Engineering,2011,9(1):83-86.

[6] 陳堅(jiān)強(qiáng),張益榮.基于Richardson插值法的CFD驗(yàn)證和確認(rèn)方法的研究 [J].空氣動(dòng)力學(xué)學(xué)報(bào),2012,30(2):176-183.

CHEN J Q,ZHANG Y R.Verification and validation in CFD based on the Richardson extrapolation method [J].Acta Aerodynamical Sinica,2012,30(2):176-183.

[7] 劉 全,王瑞利,林 忠.非嵌入多項(xiàng)式混沌方法在拉氏計(jì)算中的應(yīng)用 [J].固體力學(xué)學(xué)報(bào),2013,33(增刊):224-233.

LIU Q,WANG R L,LIN Z.Uncertainty quantification for Lagrangian computational using non-intrusive polynomial chaos [J].Chinese Journal of Solid Mechanics,2013,33(Suppl):224-233.

[8] 張 楠,沈泓萃,姚惠之.阻力和流場(chǎng)的CFD不確定度分析探討 [J].船舶力學(xué),2008,12(2):211-224.

ZHANG N,SHEN H C,YAO H Z.Uncertainty analysis in CFD for resistance and flow field [J].Journal of Ship Mechanics,2008,12(2):211-224.

[9] OBERKAMPF W,ROY C.Verification and validation in the scientific computing [M].New York:Cambridge University Press,2010.

[10] 馬智博,鄭 淼,尹建偉,等.爆轟模擬不確定度的量化方法 [J].計(jì)算物理,2011,28(1):66-74.

MA Z B,ZHENG M,YIN J W,et al.Quantification of uncertainties in detonation simulations [J].Chinese Journal of Computational Physics,2011,28(1):66-74.

[11] HELTON J,JOHNSON J,SALLABERRY C,et al.Survey of sampling-based method for uncertainty and sensitivity analysis [J].Reliab Eng Syst Safe,2006,91(11):1175-1209.

[12] HELTONl J,JOHNSON J,SALLABERRY C.Quantification of margins and uncertainties:example analyses from reactor safety and radioactive waste disposal involving the separation of aleatory and epistemic uncertainty [J].Reliab Eng Syst Safe,2008,99(1):175-195.

[13] JANSSEN H.Monte-Carlo based uncertainty analysis:sampling efficiency and sampling convergence [J].Reliab Eng Syst Safe,2013,109:123-132.

[14] 皮 霆,張?jiān)魄?吳景錸.基于多項(xiàng)式混沌方法的柔性多體系統(tǒng)不確定性分析 [J].中國(guó)機(jī)械工程,2011,22(19):2341-2348.

PI T,ZHANG Y Q,WU J L.Uncertainty analysis of flexible multibody systems using the polynomial chaos methods [J].China mechanical Engineering,2011,22(19):2341-2348.

[15] 王曉東,康 順.多項(xiàng)式混沌方法在隨機(jī)方腔流動(dòng)模擬中的應(yīng)用 [J].中國(guó)科學(xué)E輯,2011,41(6):790-798.

WANG X D,KANG S.Application of polynomial chaos on numerical simulation of stochastic cavity flow [J].Science China Series E,2011,41(6):790-798.

[16] CHEN Q,GOTTLIEB D,HESTHAVEN J.Uncertainty analysis for the steady state flows in a dual throat nozzle [J].J Comput Phys,2005,204:119-148.

[17] 水鴻壽.一維流體力學(xué)差分方法 [M].北京:國(guó)防工業(yè)出版社,1998.

SHUI H S.The difference method of one dimensional fluid mechanics [M].Beijing:National Defense Press,1998.