梁 霄,王瑞利

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

混合不確定度量化方法及其在計(jì)算流體動(dòng)力學(xué)迎風(fēng)格式中的應(yīng)用*

梁 霄1,2,王瑞利1

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

針對(duì)流體力學(xué)數(shù)值求解間斷問(wèn)題時(shí)，初始狀態(tài)含有偶然和認(rèn)知混合型的不確定性，將認(rèn)知不確定度作為外層，偶然不確定度作為內(nèi)層，分別使用非嵌入多項(xiàng)式混沌方法(non-intrusive polynomial chaos,NIPC)和概率盒(P-box)理論處理偶然不確定度和認(rèn)知不確定度，發(fā)展了流體力學(xué)數(shù)值求解過(guò)程中，初始狀態(tài)含有混合不確定度傳播量化的一種方法。以迎風(fēng)格式和黎曼解法器求解Sod問(wèn)題為例，評(píng)估了左狀態(tài)密度(偶然不確定度)和理想氣體多方指數(shù)(認(rèn)知不確定度)對(duì)模型輸出結(jié)果的影響，有效驗(yàn)證了該方法的可行性。

爆炸力學(xué);建模與模擬;非嵌入多項(xiàng)式混沌法;混合不確定度量化;概率盒

在重大復(fù)雜工程領(lǐng)域,眾多工程事故已經(jīng)引發(fā)社會(huì)各界的關(guān)注。人們總是抱著一個(gè)良好的愿望,希望能逼真建模,研制高可信度的數(shù)值模擬軟件,精確模擬各類事故的發(fā)生過(guò)程，科學(xué)分析其產(chǎn)生的機(jī)理。因此,數(shù)值模擬軟件必須能真實(shí)反映實(shí)際系統(tǒng)的外部表征、內(nèi)在特性及其動(dòng)作行為等現(xiàn)象,其數(shù)值模擬結(jié)果才能使人置信。但由于物理過(guò)程的復(fù)雜性及人們的認(rèn)知缺陷,物理建模與數(shù)值模擬始終存在誤差和不確定性。為此,要利用數(shù)值模擬結(jié)果分析機(jī)理和模擬過(guò)程,必須量化數(shù)值模擬結(jié)果的不確定性,評(píng)估建模與模擬的可信度。但長(zhǎng)期以來(lái),物理建模與數(shù)值模擬軟件的研制者重點(diǎn)關(guān)注的是軟件自身的研制和對(duì)實(shí)際問(wèn)題的順利模擬,即算下去的問(wèn)題,對(duì)軟件自身的正確性、模擬結(jié)果的不確定性和建模與模擬的可信度,即算多準(zhǔn)的問(wèn)題,一直沒(méi)有給予足夠的重視。尤其是對(duì)于理論、實(shí)驗(yàn)無(wú)法解釋或難于解決的科學(xué)問(wèn)題,或做實(shí)驗(yàn)耗資巨大的工程設(shè)計(jì),很想利用數(shù)值模擬真實(shí)、全過(guò)程,全系統(tǒng)、反復(fù)精密地進(jìn)行模擬的研究者持一種矛盾的心態(tài)：既想利用數(shù)值模擬這條快捷經(jīng)濟(jì)的重要研究途徑,又對(duì)程序的正確性和數(shù)值模擬提供的結(jié)果心存疑慮。為了給復(fù)雜裝置設(shè)計(jì)提供高效可靠的數(shù)值模擬工具,必須展開計(jì)算模型正確性和可信度的評(píng)估研究。

流體力學(xué)方程及求解方法是眾多復(fù)雜工程多物理過(guò)程建模與模擬(modeling & simulation，M&S)中的最基本模型,此計(jì)算模型的驗(yàn)證與確認(rèn)(verification and validation,V&V)受到高度重視[1-7]。W.L.Oberkampf等[8]綜述了機(jī)械工程領(lǐng)域中科學(xué)計(jì)算V&V的發(fā)展,系統(tǒng)論述了科學(xué)計(jì)算中V&V的基本概念、原理及步驟和發(fā)展過(guò)程。王瑞利等[5]綜述了V&V的相關(guān)概念、 術(shù)語(yǔ)、規(guī)范、可信度評(píng)估方法和應(yīng)用等的研究進(jìn)展,就M&S的V&V的幾個(gè)關(guān)鍵問(wèn)題進(jìn)行了概括,構(gòu)建了復(fù)雜工程M&S的V&V的知識(shí)指南。確認(rèn)是建立在定量準(zhǔn)確評(píng)價(jià)的概念上,是評(píng)價(jià)模型描述實(shí)際物理過(guò)程/實(shí)驗(yàn)(試驗(yàn))的準(zhǔn)確程度、模型形式與參數(shù)的適應(yīng)范圍,確認(rèn)刻畫實(shí)際問(wèn)題/實(shí)驗(yàn)(試驗(yàn))區(qū)域的形狀,預(yù)測(cè)可信度評(píng)估等重要活動(dòng)。對(duì)于流體力學(xué)確認(rèn)的核心是隨機(jī)輸入?yún)?shù)對(duì)模擬輸出響應(yīng)量不確定傳播進(jìn)行量化的方法。

不確定度量化(uncertainty quantification,UQ)方法是利用各種數(shù)學(xué)和統(tǒng)計(jì)學(xué)方法,對(duì)隨機(jī)因素或認(rèn)知缺陷盡可能給出分布或區(qū)間描述,建立建模和模擬中隨機(jī)或認(rèn)知不確定性因素和建模與模擬描述系統(tǒng)響應(yīng)量之間的關(guān)系,進(jìn)而量化影響系統(tǒng)性能的不確定性,給出系統(tǒng)響應(yīng)量性能分析的不確定度度量的方法。對(duì)于不確定度量化有多種方法,其中敏感度分析方法、抽樣方法、代理模型方法、隨機(jī)微分方程方法和參數(shù)率定與優(yōu)化方法已在復(fù)雜工程中得到很好的應(yīng)用。

在上述幾種方法中,隨機(jī)微分方程是20世紀(jì)中葉發(fā)展起來(lái)的一種不確定度量化方法。廣義的隨機(jī)微分方程包括隨機(jī)過(guò)程驅(qū)動(dòng)的微分系統(tǒng)和系數(shù)為隨機(jī)量的微分方程,狹義上的隨機(jī)微分方程往往專指前者。對(duì)于隨機(jī)過(guò)程驅(qū)動(dòng)的微分系統(tǒng),方程的解往往是一個(gè)隨機(jī)過(guò)程函數(shù)。雖然求解常微分方程的方法不適用于求解隨機(jī)微分方程,然而,有些復(fù)雜過(guò)程可以直接使隨機(jī)過(guò)程建模、模型的不確定性,直接表現(xiàn)在模型本身上而不是輸入?yún)?shù)的不確定性,隨機(jī)微分方程是研究不確定度量化的很好途徑。另一方面,系數(shù)為隨機(jī)量的微分方程可以直接處理為隨機(jī)參數(shù)的微分方程,常見(jiàn)的分析方法是隨機(jī)譜方法。

隨機(jī)譜方法是使用譜逼近隨機(jī)微分方程或參量,然后將其分解為獨(dú)立的確定性分量和隨機(jī)分量來(lái)量化不確定性。其數(shù)學(xué)思想是數(shù)學(xué)模型中的每一個(gè)參量利用正交多項(xiàng)式如(如Hermite多項(xiàng)式)展開成無(wú)窮級(jí)數(shù)項(xiàng),在實(shí)際應(yīng)用中取有限項(xiàng)。無(wú)窮級(jí)數(shù)第1項(xiàng)表示參量的期望,第2項(xiàng)表示高斯隨機(jī)波動(dòng),第3項(xiàng)及高階項(xiàng)表示非高斯隨機(jī)波動(dòng)。隨機(jī)譜方法已在建模和模擬不確定度量化中得到了廣泛應(yīng)用。常見(jiàn)的2種方法:KL展開(Karhunen-Loeve expansion)和多項(xiàng)式混沌方法(polynomial chaos,PC)。本文中，將著重介紹PC在Riemann問(wèn)題不確定度量化中的應(yīng)用。

1 流體力學(xué)方程組與迎風(fēng)格式

一維拉氏流體力學(xué)方程組[9]:

(1)

式中：ρ、u、E、e和p分別為流體密度、流體速度、總能、比內(nèi)能和流體壓強(qiáng),v為比容,γ為理想氣體多方指數(shù)。一維拉氏流體力學(xué)方程組迎風(fēng)型格式見(jiàn)文獻(xiàn)[9]

2 多項(xiàng)式混沌法

多項(xiàng)式混沌法其數(shù)學(xué)思想是數(shù)學(xué)模型中的每一個(gè)參量利用正交多項(xiàng)式(如Hermite多項(xiàng)式)展開成無(wú)窮級(jí)數(shù)項(xiàng)。設(shè)系統(tǒng)的解具有如下形式：

(2)

(3)

IPC是將數(shù)學(xué)模型中的每個(gè)參量進(jìn)行正交多項(xiàng)式展開,然后利用Galerkin映射,得到相應(yīng)的隨機(jī)控制方程,但是每一次展開都需要重新書寫控制方程。因此，IPC無(wú)法利用已有程序,需要對(duì)已有的數(shù)值求解程序進(jìn)行大量修改或者重新研制計(jì)算程序。

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

3 P-box的構(gòu)造方法

若問(wèn)題中同時(shí)存在偶然不確定和認(rèn)知不確定,需要將其結(jié)合起來(lái)考慮,此時(shí)可采用概率盒方法進(jìn)行描述，即將各參數(shù)分為內(nèi)外2層循環(huán)：外層循環(huán)考慮認(rèn)知不確定參數(shù),使用拉丁超立方取樣(Latin hypercube sampling,LHS)方法抽樣；內(nèi)層循環(huán)考慮偶然不確定參數(shù),取值由隨機(jī)抽樣給出。

P-box計(jì)算步驟:

(1)外層循環(huán)考慮認(rèn)知不確定參數(shù),從其區(qū)間中,采用LHS方法得到N個(gè)樣本；內(nèi)層循環(huán)考慮偶然不確定參數(shù),利用Monte Carlo(MC)抽樣得到M個(gè)樣本,作為左端密度初值;

(2)固定一個(gè)外層認(rèn)知不確定變量γ,計(jì)算出相應(yīng)的M個(gè)結(jié)果;

(3)尋找相應(yīng)物理量的最大、最小值,然后分成100份,計(jì)算每一份對(duì)應(yīng)的樣本容量,繪出相應(yīng)的累積分布曲線;

(4)對(duì)剩余的N-1個(gè)認(rèn)知不確定數(shù)施行同樣運(yùn)算,劃出N條累積分布曲線;

(5)選取最內(nèi)層和最外層的那一條作為概率邊界，構(gòu)成一個(gè)P-box;

(6)讀出概率盒的意義。

4 混合不確定量化方法在Sod問(wèn)題中的應(yīng)用

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

Sod激波管問(wèn)題是初始間斷的分解問(wèn)題,在一個(gè)無(wú)限長(zhǎng)管道中,在x=0處的左端是高壓氣體,右邊是低壓氣體,在初始時(shí)刻處于靜止?fàn)顟B(tài),兩者之間用薄膜隔離。t=0時(shí)刻薄膜破裂,形成向左稀疏波和向右激波。Sod問(wèn)題計(jì)算區(qū)域?yàn)閇-1,1],初值為:

x≤0時(shí)，ρL～N(1.0,2.5×10-7),uL=0.0,pL=1.0,γL∈[1.4,3.1];

x>0時(shí),ρR=0.125,uR=0.0,pR=0.1,γR=1.4。

采用理想氣體狀態(tài)方程:

p=(γ-1)ρe

此問(wèn)題的特點(diǎn)是左右2種狀態(tài)的γ不一樣,是一個(gè)多介質(zhì)激波管問(wèn)題,有解析解。但其中左右密度量含有不確定因素,對(duì)其量化才能實(shí)施計(jì)算模型的確認(rèn)。對(duì)于左端炸藥密度,由于加工過(guò)程中炸藥顆粒凝結(jié)的不均勻性,導(dǎo)致密度存在隨機(jī)性。

假設(shè)其服從正態(tài)分布,期望為μ=1.4 g/cm3,標(biāo)準(zhǔn)差σ=0.5 mg/cm3,即炸藥密度服從正態(tài)分布：

N(μ,σ2)=N(1.4,1.5×10-7)

圖2 多方指數(shù)分布散點(diǎn)Fig.2 Scatter of polytrophic exponent

對(duì)炸藥密度采用MC方法,取樣1 000次,分布直方圖如圖1所示。

對(duì)于多方指數(shù)γ，不知道它的任何信息,只知道它位于[1.4，3.1],在實(shí)際設(shè)計(jì)或計(jì)算中采用哪個(gè)數(shù)值由設(shè)計(jì)者決定,這種選取的過(guò)程與認(rèn)識(shí)水平有關(guān),屬于認(rèn)知不確定度。利用LHS采樣10次,如圖2所示。

當(dāng)固定γ時(shí),會(huì)得到1 000組狀態(tài)確定性Sod問(wèn)題的樣本點(diǎn),將其代入迎風(fēng)格式的流體力學(xué)程序和精確解法器,得到1 000組響應(yīng)量的樣本狀態(tài),γ=2.29時(shí),計(jì)算到t=0.22時(shí),密度、壓力、速度、能量的狀態(tài)構(gòu)成的馬尾圖,如圖3所示。從圖3可以看出,在某些地方解的分散度很明顯。當(dāng)γ變化時(shí),會(huì)得到10張類似的馬尾圖。

圖3 Sod問(wèn)題1000組計(jì)算結(jié)果Fig.3 1000 sets of computational results for the Sod problem

4.1 可信范圍

對(duì)于上述計(jì)算結(jié)果,信息量巨大,現(xiàn)在提取有用的信息。首先,將10張馬尾圖合并,得到10×1 000=10 000組響應(yīng)量樣本狀態(tài),采用NIPC方法開展不確定度量化。

圖4給出了計(jì)算得到的密度、壓力、速度、能量的期望和方差,從圖可以看出有些地方方差很大,精確求解黎曼解法器的方差比迎風(fēng)格式的高。

圖5給出了迎風(fēng)格式可信范圍和精確黎曼解法器可信范圍,從圖可以看出有80%的地方完全重疊,其他局部重疊或者不重疊,說(shuō)明采用迎風(fēng)格式求解Sod問(wèn)題是可信的。

圖4 t=0.22時(shí)Sod問(wèn)題的期望和方差Fig.4 Expectation and variance for the Sod problem at t=0.22

圖5 2種方法計(jì)算Sod可信區(qū)間的比較Fig.5 Comparison of confidential intervals from two computational methods in solving the Sod problem

4.2 Sod問(wèn)題的P-box

圖6 密度函數(shù)的累積分布函數(shù)族Fig.6 Cumulative distribution function ensemble of density function

利用第3節(jié)中的方法,構(gòu)造Sod問(wèn)題中迎風(fēng)格式的P-box,做出在x=0.1,t=0.22處,密度函數(shù)的累積分布函數(shù),如圖6所示,然后選取最內(nèi)層和最外層概率邊界曲線構(gòu)成一個(gè)P-box,如圖7所示。

從圖7可以看出，在t=0.22,x=0.1時(shí), 密度ρ(x,t)≤0.555 1的概率0.1

同時(shí),做出黎曼解法器的P-box,并與迎風(fēng)格式做比較,如圖8所示,不難看出迎風(fēng)格式構(gòu)成的P-box基本包含黎曼解法器的P-box,兩者差異不大,這說(shuō)明采用迎風(fēng)格式求解Sod問(wèn)題是有效的。

圖7 利用迎風(fēng)格式構(gòu)造的密度函數(shù)的P-boxFig.7 P-box constructed by the upwind schemefor density function

圖8 用2種方法構(gòu)成的密度函數(shù)的P-boxFig.8 P-box constructed by two methodsfor density function

5 結(jié) 論

(1)流體力學(xué)數(shù)值模擬過(guò)程中初始狀態(tài)或輸入條件常含有偶然和認(rèn)知混合型的不確定性，對(duì)模擬結(jié)果的響應(yīng)量有很大影響。

(2)針對(duì)流體力學(xué)初始狀態(tài)中含有偶然不確定性，給出了非嵌入式多項(xiàng)式混沌法量化方法。針對(duì)初始狀態(tài)中含有認(rèn)知不確定度，給出了P-box量化方法，提出了一種流體力學(xué)輸入?yún)?shù)對(duì)模擬輸出響應(yīng)量不確定性傳播量化的方法。

(3)以迎風(fēng)格式和黎曼解法器求解Sod問(wèn)題為例，給出了Sod問(wèn)題左狀態(tài)密度ρL偶然不確定度和多方指數(shù)γL認(rèn)知不確定度，分別采用非嵌入式多項(xiàng)式混沌法和P-box量化結(jié)果，給出了流體力學(xué)方程組迎風(fēng)格式數(shù)值求解可信度與輸入不確定度的范圍關(guān)系，為流體力學(xué)數(shù)值求解過(guò)程中不確定性傳播量化和可信度評(píng)估提供了一種有效方法。

[1] Roache P J. Verification of codes and calculations[J]. AIAA Journal, 1998,36(5):696-702.

[2] Staf A A. Guide for the verification and validation of computational fluid dynamics simulations[M]. Reston: American Institute of Aeronautics and Astronautics, 1998:126-180.

[3] American Nuclear Society. Guidelines for the verification and validation of scientific and engineering computer programs for the nuclear industry[M]. ANSI/ANS-10.4-1987, 1987:63-119.

[4] 鄧小剛,宗文剛,張來(lái)平,等.計(jì)算流體力學(xué)中的驗(yàn)證與確認(rèn)[J].力學(xué)進(jìn)展,2007,37(2):279-288. Deng Xiaogang, Zong Wengang, Zhang Laiping, et al. Verification and validation in computational fluid dynamics[J]. Advances in Mechanics, 2007,37(2):279-288.

[5] 王瑞利,溫萬(wàn)治.復(fù)雜工程建模與模擬的驗(yàn)證與確認(rèn)[J].計(jì)算機(jī)輔助工程,2014,23(4):61-68. Wang Ruili, Wen Wanzhi. Advances in verification and validation of modeling and simulation of the complex engineering[J]. Computer Aided Engineering, 2014,23(4):61-68.

[6] 劉全,王瑞利,林忠,等.爆轟計(jì)算JWL狀態(tài)方程參數(shù)不確定度研究[J].爆炸與沖擊,2013,33(6):647-654. Liu Quan, Wang Ruili, Lin Zhong, et al. Uncertainty quantification for JWL EOS parameters in explosive numerical simulation[J]. Explosion and Shock Waves, 2013,33(6):647-654.

[7] 劉全,王瑞利,林忠.非嵌入式多項(xiàng)式混沌方法在拉氏計(jì)算中的應(yīng)用[J].固體力學(xué)學(xué)報(bào), 2013,33(增):224-233. Liu Quan, Wang Ruili, Lin Zhong. Uncertainty quantification for Lagrange computation using non-intrusive polynomial chaos[J]. Chinese Journal of Solid Mechanics, 2013,33(suppl):224-233.

[8] Oberkampf W L, Roy C. Verification and validation in scientific computing[M]. New York: Cambridge University Press, 2010:241-305.

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

[10] Xiu Dongbin, Karniadkis G E. The Wiener-Askey polynomial chaos for stochastic differential equations[J]. SIAM Journal on Scientific Computing, 2002,24(2):619-644.

[11] Cameron R H, Martin W T. The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals[J]. Annals of Mathematics, 1947,48(2):385-398.

(責(zé)任編輯 張凌云)

Mixed uncertainty quantification and its application in upwind scheme for computational fluid dynamics (CFD)

Liang Xiao1,2, Wang Ruili1

(1.InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094,China;2.CollegeofMathematics,ShandongUniversityofScienceandTechnology,Qingdao266590,Shandong,China)

Both aleatory uncertainty and epistemic uncertainty exist in the initial and boundary conditions when we numerically solve the CFD with sharp discontinuity. In this paper, mixed uncertainty quantification approaches are developed to deal with this situation. Specifically, the outer level uncertainty is linked to epistemic uncertainties, and the inner uncertainty is linked to aleatory uncertainty. The non-intrusive polynomial chaos method is utilized to cope with the aleatory uncertainties, while the P-box theory is used to deal with the epistemic uncertainties, and the upwind scheme and Riemann solver are used to solve the deterministic system. We apply this method to the Sod problem in the CFD, and acquire preferable effect. This method evaluates the influence of input uncertainty such as density (aleatory uncertainty) and polytrophic exponent (epistemic uncertainty) on the output uncertainty, the efficiency of this method is also proved. This method is also helpful in evaluating the degree of confidence and validation of the result from modeling and simulation by other models.

mechanics of explosion; modeling and simulation; non-intrusive polynomial chaos; mixed uncertainty quantification; P-box

10.11883/1001-1455(2016)04-0509-07

2014-12-24；

2015-03-27

國(guó)家自然科學(xué)基金項(xiàng)目(11372051,11475029);中國(guó)工程物理研究院科學(xué)基金項(xiàng)目(2015B0202045);山東科技大學(xué)人才引進(jìn)項(xiàng)目(2013RCJJ027)

梁 霄(1984— )，男，博士，講師，mathlx@163.com。

O385國(guó)標(biāo)學(xué)科代碼：13035

A