流體力學仿真軟件可信度評估與預測能力

王瑞利++梁霄++喻強

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

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

中圖分類號： TB115.1文獻標志碼： 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

基金項目： 國家自然科學基金（11372051，91630312，11475029）；國防基礎科研計劃（C1520110002）；

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

作者簡介： 王瑞利（1964—），男，研究員，研究方向為計算流體力學及其應用軟件，（Email）wang_ruili@iapcm.ac.cn0引言

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

1爆轟流體動力學模型及其不確定性

1.1爆轟流體力學物理模型

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

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

圓筒試驗結構示意和計算模型見圖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點起爆或者面爆，炸藥采用JWL狀態(tài)方程和Wilkins反應率模型，紫銅采用Gruneisen狀態(tài)方程和SG本構模型。對每個因素采用抽樣技術，然后將其輸入到自主研發(fā)的確定性軟件LAD2D[11]中，產生敏感度分析響應量的樣本。各因素樣本計算值見表2。計算參數(shù)除因素欄說明外，其余均按統(tǒng)一計算條件取值：JWL參數(shù)取R1=4.6，R2=1.3，ω=0.38；燃燒函數(shù)中的參數(shù)取nb=1.3，rb=2.1，采用壓縮比起爆時取σ=1.03，采用時間起爆時按惠更斯原理預先計算起爆時間。

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

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

由于JWL狀態(tài)方程中參數(shù)A，B，C依賴于參數(shù)R1，R2，w，為此，表征或標定參數(shù)R1，R2，w非常重要。對于大多數(shù)炸藥，憑先驗估計這3個參數(shù)屬于認知不確定度，R1∈[4.0，5.0]，R2∈[1.0，2.0]，w∈[0.2，0.4]?；诶》匠闃臃椒?，本文發(fā)展高維抽樣技術，即ξi=a+Δx·θ，對R1，R2，w這3個不確定性參數(shù)進行抽樣，其中100組樣本和1 000組樣本抽樣結果見圖3。

考慮到爆轟產物的壓力恒為正，則有不等式約束條件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]，可計算出A，B和C這3個參數(shù)，再通過A>0，B>0，C>0，就能判斷樣本的有效性。參數(shù)R1，R2和w在100組樣本點和1 000組樣本點中的有效樣本見圖4。在100組樣本中有81組樣本有效，在1 000組樣本中有840組樣本有效。有效樣本求出的參數(shù)A，B和C見圖5。圖中為A>0，B>0，C>0的有效樣本，滿足物理要求。

3.2響應量徑向壁位置和速度

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

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

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

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

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

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

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

3.4Wilkins反應率不確定性量化及可信度評估

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

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

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

4結束語

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

本文只是嘗試模型確認的思路，通過實驗說明其可行性。在復雜工程建模與模擬中，不確定性涉及內容較多，需要開展大量研究工作，應引起工程仿真學者的高度重視。

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

參考文獻：

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

[2]王瑞利， 梁霄. 爆轟數(shù)值模擬中物理模型分層確認實驗研究[J]. 中國測試， 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ù)學方法[J]. 中國科學： 數(shù)學， 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]王瑞利， 溫萬治. 復雜工程建模和模擬的驗證與確認[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]. 力學進展， 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]王瑞利， 劉全， 溫萬治. 非嵌入式多項式混沌法在爆轟產物JWL參數(shù)評估中的應用[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]梁霄， 王瑞利. 爆轟流體力學模型敏感度分析與模型確認[J]. 物理學報， 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]王瑞利， 林忠， 溫萬治， 等. 多介質拉氏自適應流體動力學軟件LAD2D研制及其應用[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做功能力的實驗研究[J]. 高壓物理學報， 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]于川， 劉文翰， 李良忠， 等. 鈍感炸藥圓筒試驗與爆轟產物JWL狀態(tài)方程研究[J]. 高壓物理學報， 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]顏王吉， 曹詩澤， 任偉新. 結構系統(tǒng)識別不確定性分析的Bayes方法及其進展[J]. 應用數(shù)學和力學， 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]何炎高， 徐定華， 陳瑞林. 紡織材料設計反問題的貝葉斯統(tǒng)計推斷方法[J]. 紡織學報， 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月計 算 機 輔 助 工 程Computer Aided EngineeringVol.26 No.6Dec. 2017endprint