流體力學(xué)仿真軟件可信度評(píng)估與預(yù)測(cè)能力

王瑞利++梁霄++喻強(qiáng)

摘要： 針對(duì)爆轟流體力學(xué)過程的控制方程和輔助方程，剖析模型參數(shù)、模型形式和數(shù)值方法等的不確定因素，列舉爆轟流體力學(xué)模型的不確定性因素并開展敏感度分析。針對(duì)幾個(gè)重要參數(shù)，發(fā)展高維參數(shù)抽樣技術(shù)，結(jié)合爆轟產(chǎn)物JWL狀態(tài)方程的圓筒試驗(yàn)，開展自主研發(fā)的爆轟彈塑性流體動(dòng)力學(xué)LAD2D仿真軟件的可信度評(píng)估研究，給出可信度評(píng)估結(jié)果。此方法和研究思路為復(fù)雜工程仿真軟件可信度評(píng)估提供一種有效手段。

關(guān)鍵詞： 爆轟過程； 唯象模型； 數(shù)值模擬； 不確定性量化； 敏感度分析； 高維抽樣

中圖分類號(hào)： TB115.1文獻(xiàn)標(biāo)志碼： B

Confidence evaluation and prediction ability of

hydrodynamics simulation software

WANG Ruili1， LIANG Xiao2， YU Qiang3

（1. Beijing Institute of Applied Physics and Computational Mathematics， Beijing 100094， China；

2. College of Mathematics and Systems Science， Shandong University of Science and Technology， Qingdao 266590，

Shandong， China； 3. Beijing Anwise Technology Co.， Ltd.， Beijing 100022， China）

Abstract： As to the control equations and auxiliary equations of detonation fluid dynamics， the uncertain factors， such as the model parameters， the model form， and numerical methods， are analyzed， the uncertainty factors of detonation fluid dynamics model are listed， and their sensitivities are analyzed. Focusing on several important parameters， the high dimension parameter sampling method is developed. The reliability evaluation of the LAD2D software for the simulation of detonation fluid dynamics is carried out by combing with the JWL equation of state in the cylinder test. The reliability evaluation results are given finally. The method and research idea provides a very effective technology for evaluating the credibility of complex engineering simulation software.

Key words： detonation process； phenomenological model； numerical simulation； uncertainty quantification； sensitivity analysis； high dimension sampling

收稿日期： 2017[KG*9〗08[KG*9〗24修回日期： 2017[KG*9〗08[KG*9〗29

基金項(xiàng)目： 國家自然科學(xué)基金（11372051，91630312，11475029）；國防基礎(chǔ)科研計(jì)劃（C1520110002）；

中國工程物理研究院科學(xué)基金（2015B0202045）

作者簡(jiǎn)介： 王瑞利（1964—），男，研究員，研究方向?yàn)橛?jì)算流體力學(xué)及其應(yīng)用軟件，（Email）wang_ruili@iapcm.ac.cn0引言

高能炸藥爆轟和可燃?xì)怏w燃燒爆炸過程都是多尺度、多物理場(chǎng)耦合的復(fù)雜系統(tǒng)，涉及高溫、高壓、高速以及材料相變等多種介質(zhì)相互作用、相互混合，使理論和實(shí)驗(yàn)研究困難增大。[12]為此，國家高端領(lǐng)域復(fù)雜系統(tǒng)可靠性認(rèn)證、大型裝置性能評(píng)估和民用設(shè)備意外爆炸事故分析與預(yù)防等成為科學(xué)難題和瓶頸。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，高置信度建模、高可信度應(yīng)用軟件開發(fā)與數(shù)值仿真技術(shù)逐漸彌補(bǔ)理論和實(shí)驗(yàn)的缺陷，成為理論、實(shí)驗(yàn)、模擬與科學(xué)、技術(shù)和工程相關(guān)領(lǐng)域交叉研究的重要課題。復(fù)雜系統(tǒng)的數(shù)值仿真技術(shù)是信息時(shí)代世界各國特別是發(fā)達(dá)國家激烈競(jìng)爭(zhēng)的技術(shù)制高點(diǎn)，是一個(gè)國家綜合國力、科技創(chuàng)新力和國防裝備戰(zhàn)斗力的重要標(biāo)志，是未來提高武器裝備研制的有效手段。但是，高能炸藥爆轟和可燃?xì)怏w燃燒爆炸過程具有瞬態(tài)性，其在極短時(shí)間內(nèi)劇烈變化且各種因素相互耦合，難以獨(dú)立（分解）處理，理論和實(shí)驗(yàn)難度都很高。基于對(duì)此類問題認(rèn)知的缺陷，很多研究只能針對(duì)所關(guān)心的重大核心問題，忽略次要因素，通過適當(dāng)?shù)暮?jiǎn)化處理，建立適于解決某一問題的物理模型或唯象模型，探尋不斷完善、逐漸逼近實(shí)際的數(shù)值仿真技術(shù)。這些研究方法最大的缺陷是由于模型參數(shù)、模型形式、逼近方法等眾多不確定性融合在一起形成的強(qiáng)不確定性或模擬預(yù)測(cè)結(jié)果的隨機(jī)性。[34]因此，國防安全領(lǐng)域中先進(jìn)新型武器設(shè)計(jì)、現(xiàn)役武器維護(hù)、退役武器處理、可燃?xì)怏w爆炸等數(shù)值仿真技術(shù)都面臨強(qiáng)不確定性或隨機(jī)性的挑戰(zhàn)，現(xiàn)有方法和技術(shù)還不能處理此類強(qiáng)不確定性問題，亟待有所突破與創(chuàng)新。近幾年，美國從國家層面推出ASC等系列計(jì)劃，其目的是推動(dòng)數(shù)字化設(shè)計(jì)和計(jì)算能力的提升，應(yīng)對(duì)來自其他國家的競(jìng)爭(zhēng)壓力。驗(yàn)證與確認(rèn)已成為提高可信度的最佳途徑，其最大瓶頸仍是模型參數(shù)、模型形式、逼近方法等眾多因素融合在一起的強(qiáng)不確定性，嚴(yán)重影響仿真軟件的可信度。[57]中國亟需從戰(zhàn)略高度認(rèn)識(shí)到國防或CAE領(lǐng)域數(shù)值仿真技術(shù)中強(qiáng)不確定性量化和可信度評(píng)估的重要性。實(shí)際上，不論在實(shí)際生活還是科學(xué)研究中，研究人員都會(huì)面臨對(duì)某些系統(tǒng)的運(yùn)行性能進(jìn)行評(píng)價(jià)，得到的評(píng)價(jià)結(jié)果對(duì)妥善使用系統(tǒng)或者改進(jìn)系統(tǒng)至關(guān)重要。endprint

1爆轟流體動(dòng)力學(xué)模型及其不確定性

1.1爆轟流體力學(xué)物理模型

爆轟流體力學(xué)物理模型是由雙曲型的流體力學(xué)（偏微分方程組）與炸藥唯象反應(yīng)率模型（一階常微分方程）、物態(tài)方程（復(fù)雜函數(shù)關(guān)系式）耦合在一起的非線性偏微分方程組。

2.2基于圓筒試驗(yàn)的爆轟模型多因素敏感度分析

圓筒試驗(yàn)結(jié)構(gòu)示意和計(jì)算模型見圖1。炸藥取TNT炸藥，性能參數(shù)為γT=3.1，ρT=1.634 g/cm3，DCJ=6.932 km/s；紫銅物理性能參數(shù)為γC=3.68，ρC=8.93 g/cm3，c=3.94 km/s。在O點(diǎn)起爆或者面爆，炸藥采用JWL狀態(tài)方程和Wilkins反應(yīng)率模型，紫銅采用Gruneisen狀態(tài)方程和SG本構(gòu)模型。對(duì)每個(gè)因素采用抽樣技術(shù)，然后將其輸入到自主研發(fā)的確定性軟件LAD2D[11]中，產(chǎn)生敏感度分析響應(yīng)量的樣本。各因素樣本計(jì)算值見表2。計(jì)算參數(shù)除因素欄說明外，其余均按統(tǒng)一計(jì)算條件取值：JWL參數(shù)取R1=4.6，R2=1.3，ω=0.38；燃燒函數(shù)中的參數(shù)取nb=1.3，rb=2.1，采用壓縮比起爆時(shí)取σ=1.03，采用時(shí)間起爆時(shí)按惠更斯原理預(yù)先計(jì)算起爆時(shí)間。

基于上述樣本數(shù)據(jù)，爆轟壓力、爆轟速度和爆轟位置的敏感性分析結(jié)果見圖2。各因素對(duì)爆轟壓力、爆轟位置影響的敏感性差異不大，而對(duì)爆轟速度敏感性差別比較大。將12個(gè)因素按速度敏感性進(jìn)行排序，見表3。從圖2和表3可以看出，對(duì)于爆轟模型，各影響因素的敏感性大小排序?yàn)椋簳r(shí)間起爆燃燒下γb，體積起爆燃燒下γb，時(shí)間起爆下JWL ω，體積起爆燃燒下nb，時(shí)間起爆下JWL R1，體積起爆下JWL ω，體積起爆下JWL R2，時(shí)間起爆下JWL R2，體積起爆下JWL R1，時(shí)間起爆燃燒下nb，體積起爆閾值σ，起爆方式。由此可知，燃燒函數(shù)和JWL產(chǎn)物狀態(tài)方程是2個(gè)敏感性強(qiáng)的模型。為此，需要重視其參數(shù)選取及模型形式。a） 爆轟壓力3產(chǎn)物狀態(tài)方程不確定性量化及軟件可信度評(píng)估通過上述敏感性分析可以看出在爆轟物理模型和數(shù)值模擬中哪些因素“很敏感”或“很關(guān)鍵”，針對(duì)“很敏感”因素開展不確定性量化，以評(píng)估仿真軟件的可信度。這里選取炸藥爆轟產(chǎn)物JWL狀態(tài)方程的唯象模型形式中的不確定參數(shù)為不確定度量化對(duì)象。

3.1產(chǎn)物狀態(tài)方程參數(shù)區(qū)間抽樣產(chǎn)生有效樣本

由于JWL狀態(tài)方程中參數(shù)A，B，C依賴于參數(shù)R1，R2，w，為此，表征或標(biāo)定參數(shù)R1，R2，w非常重要。對(duì)于大多數(shù)炸藥，憑先驗(yàn)估計(jì)這3個(gè)參數(shù)屬于認(rèn)知不確定度，R1∈[4.0，5.0]，R2∈[1.0，2.0]，w∈[0.2，0.4]。基于拉丁方抽樣方法，本文發(fā)展高維抽樣技術(shù)，即ξi=a+Δx·θ，對(duì)R1，R2，w這3個(gè)不確定性參數(shù)進(jìn)行抽樣，其中100組樣本和1 000組樣本抽樣結(jié)果見圖3。

考慮到爆轟產(chǎn)物的壓力恒為正，則有不等式約束條件A>0且B>0且C>0。根據(jù)此約束條件，以JOB9003炸藥為例[12]，從物理層面上判斷樣本的有效性。JOB9003炸藥的物理性能參數(shù)為ρ0=1.849 g/cm3，k=2.99，DCJ=8.712 km/s。由A，B和C參數(shù)與R1，R2和w參數(shù)的線性方程組[3]，可計(jì)算出A，B和C這3個(gè)參數(shù)，再通過A>0，B>0，C>0，就能判斷樣本的有效性。參數(shù)R1，R2和w在100組樣本點(diǎn)和1 000組樣本點(diǎn)中的有效樣本見圖4。在100組樣本中有81組樣本有效，在1 000組樣本中有840組樣本有效。有效樣本求出的參數(shù)A，B和C見圖5。圖中為A>0，B>0，C>0的有效樣本，滿足物理要求。

3.2響應(yīng)量徑向壁位置和速度

在眾多爆轟產(chǎn)物狀態(tài)方程的形式中，JWL狀態(tài)方程是一種不顯含化學(xué)反應(yīng)、由試驗(yàn)方法確定參數(shù)的半經(jīng)驗(yàn)狀態(tài)方程，能比較精確地描述爆轟產(chǎn)物的膨脹驅(qū)動(dòng)做功過程。圓筒試驗(yàn)是指將炸藥放入等壁厚的銅質(zhì)圓筒中，從圓筒的一端將其引爆，利用高速轉(zhuǎn)鏡式掃描相機(jī)記錄筒壁在爆轟產(chǎn)物驅(qū)動(dòng)下的膨脹過程，是確定炸藥爆轟產(chǎn)物JWL狀態(tài)方程和評(píng)估炸藥做功能力的標(biāo)準(zhǔn)化試驗(yàn)，國內(nèi)外廣泛應(yīng)用。圓筒試驗(yàn)布局見圖6，采用狹縫掃描和VISAR聯(lián)合測(cè)試方法。主炸藥以JOB9003為例，直徑為25.4 mm，長(zhǎng)度為305.0 mm；圓筒內(nèi)徑為25.4 mm，壁厚為2.6 mm，材料為紫銅。狹縫掃描采用平行光后照明技術(shù)，狹縫位置距離起爆端200.0 mm。相機(jī)轉(zhuǎn)速為6×104 r/min，對(duì)應(yīng)的掃描速度為3 km/s。VISAR測(cè)試與狹縫掃描光路垂直，光纖探頭前段距裝置外表面75.0 mm，光路與靜態(tài)裝置外表面垂直。VISAR測(cè)試位置與狹縫掃描位置對(duì)應(yīng)，距起爆端200.0 mm。同時(shí)，在圓筒末端的定常滑移爆轟段采用電探針法測(cè)定炸藥爆速。[1213]

在參數(shù)R1，R2和w的樣本點(diǎn)上進(jìn)行確定性數(shù)值計(jì)算時(shí)采用LAD2D軟件。[11]LAD2D是北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所自主開發(fā)的大型爆轟彈塑性流體動(dòng)力學(xué)程序，空間離散采用拉氏有限體積格式，能夠準(zhǔn)確捕捉自由面速度，時(shí)間離散采用預(yù)報(bào)校正方法。計(jì)算過程采用能夠有效處理拉氏計(jì)算過程中網(wǎng)格大變形問題的網(wǎng)格鄰域可變技術(shù)。為有效模擬激波傳播及其相互作用的問題，采用人工黏性方法，綜合考慮流體彈塑性和炸藥反應(yīng)過程。大量數(shù)值模擬計(jì)算結(jié)果表明，LAD2D具有較為完整的解決工程實(shí)際問題的模擬能力。在81組有效樣本上模擬計(jì)算圓筒試驗(yàn)，得到離起爆端200.0 mm狹縫面的徑向壁位置和速度曲線，見圖7。

a） 徑向壁位置b） 徑向壁速度圖 7LAD2D計(jì)算得到的基于有效參數(shù)樣本的響應(yīng)量

3.3JWL狀態(tài)方程不確定性量化及可信度評(píng)估

基于上述響應(yīng)量的樣本，通過統(tǒng)計(jì)方法，得到位置與速度隨時(shí)間演變的期望與方差。參數(shù)R1，R2和w的不確定性對(duì)圓筒試驗(yàn)離起爆端200.0 mm狹縫處徑向壁位置和速度傳播的影響見圖8。從圖7可以看出，響應(yīng)量（徑向壁位置和徑向壁速度）分散度比較大。采用多項(xiàng)式混沌方法給出的模擬結(jié)果的不確定性見圖8，將其與確認(rèn)試驗(yàn)（含有不確定性評(píng)估）的測(cè)試數(shù)據(jù)對(duì)比，就可以給出模型和軟件的可信度。同時(shí)，采用貝葉斯原理[1415]，可反推給出爆轟唯象模型的可信參數(shù)，達(dá)到高置信度的預(yù)測(cè)能力。endprint

a） 徑向壁位置b） 徑向壁速度圖 8基于有效樣本的響應(yīng)量不確定性量化

Fig.8Uncertainty quantification of response quantity based on effective parameter samples

3.4Wilkins反應(yīng)率不確定性量化及可信度評(píng)估

針對(duì)Wilkins反應(yīng)率唯象模型中的參數(shù)nb和rb，憑先驗(yàn)估計(jì)，這2個(gè)參數(shù)屬于認(rèn)知不確定度，nb∈[1.0，6.0]，γb∈[2.0，3.0]，基于拉丁方抽樣方法給出30個(gè)樣本。對(duì)應(yīng)的狹縫徑向壁位置和速度的模擬結(jié)果見圖9，采用統(tǒng)計(jì)分析方法給出模擬結(jié)果的不確定性量化見圖10。

a） 徑向壁位置b） 徑向壁速度圖 9LAD2D計(jì)算得到的參數(shù)樣本的響應(yīng)量

從圖9和10可以看出，Wilkins反應(yīng)率唯象模型分散度比JWL狀態(tài)方程小，說明此模型不確定性小。

4結(jié)束語

針對(duì)炸藥爆轟唯象模型（爆轟產(chǎn)物JWL狀態(tài)方程）中的不確定性參數(shù)，通過拉丁方抽樣技術(shù)，并結(jié)合爆轟流體力學(xué)軟件LAD2D開展圓筒試驗(yàn)計(jì)算，建立不確定性參數(shù)與響應(yīng)量的樣本，通過貝葉斯理論，校準(zhǔn)與確認(rèn)計(jì)算模型的參數(shù)。

本文只是嘗試模型確認(rèn)的思路，通過實(shí)驗(yàn)說明其可行性。在復(fù)雜工程建模與模擬中，不確定性涉及內(nèi)容較多，需要開展大量研究工作，應(yīng)引起工程仿真學(xué)者的高度重視。

針對(duì)國防或CAE領(lǐng)域機(jī)理建模與模擬中強(qiáng)不確定性現(xiàn)象密切相關(guān)的科學(xué)規(guī)律、數(shù)學(xué)方法與理論，發(fā)展復(fù)雜工程仿真軟件可信度評(píng)估方法及應(yīng)用的前沿核心技術(shù)，突破數(shù)值仿真技術(shù)未考慮不確定性的桎梏，建立國防或CAE領(lǐng)域復(fù)雜系統(tǒng)數(shù)值模擬技術(shù)中強(qiáng)不確定性高效量化方法及分析軟件，發(fā)展一批強(qiáng)預(yù)測(cè)能力的仿真軟件，是未來國防或CAE領(lǐng)域發(fā)展的重要方向。

參考文獻(xiàn)：

[1]孫錦山， 朱建士. 理論爆轟物理[M]. 北京： 國防工業(yè)出版社， 1995.

[2]王瑞利， 梁霄. 爆轟數(shù)值模擬中物理模型分層確認(rèn)實(shí)驗(yàn)研究[J]. 中國測(cè)試， 2016， 42（10）： 1320.

WANG R L， LIANG X. Research on validation experiment hierarchy of validation for physical modeling in numerical simulation of detonation[J]. China Measuement & Testing Technology， 2016， 42（10）： 1320.

[3]王瑞利， 江松. 多物理耦合非線性偏微分方程與數(shù)值解不確定度量化數(shù)學(xué)方法[J]. 中國科學(xué)： 數(shù)學(xué)， 2015， 45（6）： 723738. DOI： 10.1360/N01201400115.

WANG R L， JIANG S. Mathematical methods for uncertainty quantification in nonlinear multiphysics systems and their numerical simulations[J]. Science China： Math， 2015， 45（6）： 723738. DOI： 10.1360/N01201400115.

[4]OBERKAMPF W L， ROY C J. Verification and validation in scientific computing[M]. Cambridge： Cambridge University Press， 2010.

[5]王瑞利， 溫萬治. 復(fù)雜工程建模和模擬的驗(yàn)證與確認(rèn)[J]. 計(jì)算機(jī)輔助工程， 2014， 23（4）： 6168. DOI： 10.13340/j.cae.2014.03.013.

WANG R L， WEN W Z， Advances in verification and validation of modeling and simulation of the complex engineering[J]. Computer Aided Engineering， 2014， 23（4）： 6168. DOI： 10.13340/j.cae.2014.03.013.

[6]鄧小剛， 宗文剛， 張來平， 等. 計(jì)算流體力學(xué)中的驗(yàn)證與確認(rèn)[J]. 力學(xué)進(jìn)展， 2007， 37（2）： 279288.

DENG X G， ZONG W G， ZHANG L P， et al. Verification and validation in computational fluid dynamics[J]. Advances in Mechanics， 2007， 37（2）： 279288.

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

WANG R L， LIU Q， WEN W Z. Nonintrusive polynomial chaos methods and its application in the parameters assessment of explosion product JWL[J]. Explosion and Shock Waves， 2015， 35（1）： 915.

[8]梁霄， 王瑞利. 爆轟流體力學(xué)模型敏感度分析與模型確認(rèn)[J]. 物理學(xué)報(bào)， 2017， 66（11）： 116401. DOI： 10.7498/aps.66.116401.endprint

LIANG X， WANG R L. Sensitivity analysis and validation of detonation computational fluid dynamics model[J]. Acta Physica Sinica， 2017， 66（11）： 116401. DOI： 10.7498/aps.66.116401.

[9]HELTON J C. Conceptual and computational basis for the quantification of margins and uncertainty： SAND20093055[R].

[10]KARNIADAKIS G E， GLIMM J. Preface： uncertainty quantification in simulation science[J]. Journal of Computational Physics， 2006， 217（1）： 14.

[11]王瑞利， 林忠， 溫萬治， 等. 多介質(zhì)拉氏自適應(yīng)流體動(dòng)力學(xué)軟件LAD2D研制及其應(yīng)用[J]. 計(jì)算機(jī)輔助工程， 2014， 23（2）： 17.DOI： 10.13340/j.cae.2014.02.001.

WANG R L， LIN Z， WEN W Z， et al. Development and application of adaptive multimedia Lagrangian fluid dynamics software LAD2D[J]. Computer Aided Engineering， 2014， 23（2）： 17.DOI： 10.13340/j.cae.2014.02.001.

[12]徐輝， 孫占峰. 鈍感高能炸藥JB9014做功能力的實(shí)驗(yàn)研究[J]. 高壓物理學(xué)報(bào)， 2013， 27（4）： 582586.

XU H， SUN Z F. An experimental study on the capacity for work of insensitive high explosive[J]. Chinese Journal of High Pressure Physics， 2013， 27（4）： 582586.

[13]于川， 劉文翰， 李良忠， 等. 鈍感炸藥圓筒試驗(yàn)與爆轟產(chǎn)物JWL狀態(tài)方程研究[J]. 高壓物理學(xué)報(bào)， 1997， 11（3）： 227233.

YU C， LIU W H， LI L Z， et al. Studies on cylinder test and JWL equation of state of detonation product for insensitive high explosive[J]. Chinese Journal of High Pressure Physics， 1997， 11（3）： 227233.

[14]顏王吉， 曹詩澤， 任偉新. 結(jié)構(gòu)系統(tǒng)識(shí)別不確定性分析的Bayes方法及其進(jìn)展[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2017， 38（1）： 4459.

YAN W J， CAO S Z， REN W X. Uncertainty quantification for system identification utilizing the Bayes theory and its recent advances[J]. Applied Mathematics and Mechanics， 2017， 38（1）： 4459.

[15]何炎高， 徐定華， 陳瑞林. 紡織材料設(shè)計(jì)反問題的貝葉斯統(tǒng)計(jì)推斷方法[J]. 紡織學(xué)報(bào)， 2015， 36（1）： 2329.

HE Y G， XU D H， CHEN R L. Bayesian statistical inference method for inverse problems of textile material design[J]. Journal of Textile Research， 2015， 36（1）： 2329.（編輯武曉英）第26卷 第6期2017年12月計(jì) 算 機(jī) 輔 助 工 程Computer Aided EngineeringVol.26 No.6Dec. 2017endprint