基于多項式混沌方法對C-J 爆轟參數不確定度的分析*

梁 霄，王瑞利，胡星志，陳江濤

（1. 山東科技大學數學與系統科學學院，山東 青島 266590；2. 北京應用物理與計算數學研究所，北京 100094；3. 中國空氣動力研究與發展中心，四川 綿陽 621000）

爆轟是極其劇烈的物理化學反應，炸藥爆炸過程擴展的速度高達每秒數千米到萬米之間，所形成的溫度約為3 000～5 000 ℃，壓力高達102～104MPa，發生在極短的時間（10-9s 量級）和極窄的空間（10-4m量級），釋放極高的能量（102GW/cm2量級）[1-5]。爆轟過程的極端特征給試驗精確測量和理論表征帶來了極大挑戰。由于爆轟試驗成本高和化學反應的不穩定性，只有經過多次試驗才能標定感興趣量（quantity of interest, QoI）的取值。即使是常規炸藥，也只有有限的試驗數據。與試驗相比，建模與模擬（modeling and simulation, M&S）是一種經濟、安全的途徑。為降低運算成本，提高計算可行性，通常還使用假設、簡化和近似等手段，以降低模型的復雜性[1,4-6]。Chapman-Jouguet（C-J）理論假設流動是一維的，不考慮熱傳導、熱輻射以及黏滯摩擦等耗散效應，忽略沖擊起爆過渡和不定常效應，將復雜非線性多物理爆轟過程簡化為一個具有沖擊壓縮間斷面的一維問題。C-J 理論能建立波后QoIs 和爆速、初始密度以及唯象參數的函數關系式，能預測爆壓和爆熱等不易測量或測量成本較高的物理量，具有實用性強、適用范圍廣的優點，是高能炸藥爆轟研究的有效工具[5-6]。

C-J 理論是學術界公認的成功理論，其模擬過程中使用的物理量和參數一直被認為是確定的。事實上，由于炸藥組分的異質性以及物理量測量的隨機特征，物理量和唯象參數都會受到不確定性的干擾。例如，爆速是描述炸藥力學行為的基本物理量，某些HMX 基高能炸藥的爆速超過9×103m/s，速度的瞬時極大變化率導致缺乏精密的試驗裝置進行測試。傳統的試驗裝置如圓筒試驗，所用金屬的純度和密度不均勻性、加工公差的存在以及試驗數據的擬合誤差，都會導致測量結果的不確定性。同時，含能材料組分的細微變化、幾何邊界的不精確、極端荷載下金屬延展性變化以及黏結劑的存在也會導致試驗結果的波動[7-9]。

密度是決定爆轟性能的另一個基本物理量。密度的不確定度來自于炸藥加工過程中產生的空腔、孔隙、裂紋、扭結和顆粒結晶的隨機性以及雜質的混入[10]。此外，存儲溫度的變化也會引起密度的波動[11]。密度的微小變化會激發感度和熱點的變化，甚至會引起爆轟系統輸出結果劇烈的、奇異的甚至目前理論無法解釋的變動[7,12]。

綜上所述，爆轟實際過程中，如爆速和初始密度等物理量都會受到噪聲擾動。這類不確定度是炸藥的固有屬性，即使增加樣本容量，也很難減少或降低這種不確定性。另一方面，爆轟產物狀態方程、反應率函數等爆轟唯象模型中的參數，無法通過試驗直接標定，需要憑借工程技術人員的經驗[3,13-14]，這類不確定性也會嚴重影響炸藥爆轟數值模擬的精度和可信性。

事實上，工程技術人員和學者們都已經意識到爆轟中的不確定因素[1,5-6,9,15]。文獻[1, 6]甚至將物理量的容許不確定度（誤差）定為 ±5 %。然而，很少有學者專門研究不確定因素對QoI 的影響。實際上，不確定因素的存在使得研究人員一直面臨著試驗成本太高和數值模擬不穩定的矛盾心態。導致研究人員既想利用C-J 理論這種經濟、安全的手段，又對C-J 理論的使用范圍和預測結果心存疑慮。因此，爆轟不確定度量化（uncertainty quantification, UQ）研究是增強炸藥爆轟唯象模型的可靠性、標定模型的使用范圍、緩和試驗與數值模擬矛盾的重要途徑，具有重要的學術意義和實用價值。

近10 年來，歐美UQ 技術在水文、地理、預報、經濟、自動控制和結構力學分析等領域蓬勃發展。但爆轟UQ 研究結果較少，面臨眾多挑戰。首先，在刻畫炸藥爆轟模型中輸入不確定度方面，數理統計理論完備，技術成熟，是量化不確定度的最自然的工具。根據大數定律和中心極限定理，正態分布是描述輸入不確定度的有效工具，且參數易于通過經典假設檢驗法標定。然而，正態分布無法保證密度和爆速的非負性。其次，在如何準確得到炸藥產物狀態方程（equation of state, EOS）和反應率函數等唯象模型中參數的可信取值范圍方面，通常根據工程技術人員的經驗給定取值范圍，存在很大的人為因素。若假設參數服從此區間上的、無信息的（un-informative）均勻分布，很多重要信息會丟失。獨立同分布是數理統計中的重要假設，但爆轟系統的輸入不確定度未必滿足此假設。除此之外，針對如何評估輸入不確定度對系統響應量的影響方面，即使通過大規模工程計算，也很難給出QoIs 的概率密度函數（probability density function, PDF）的顯式函數關系。幸運的是，計算技術的快速發展和數值分析方法的進步提供了計算PDF 的有效方法。與Monte Carlo 方法相比，多項式混沌（polynomial chaos, PC）收斂速度快、運算成本低，是目前工程技術領域處理不確定度傳播的主流工具[3,16-18]。特別是非嵌入式PC，不需要更改程序，且能給出Gauss 統計量的顯式表達，是處理爆轟UQ 的潛在工具。但是在爆轟UQ 研究中，會面臨輸入不確定度概率類型不一致、隨機變量不獨立的問題，導致PC 理論無法直接應用。

面對眾多挑戰，本文中，結合真實的試驗數據，通過參數估計和Anderson-Darling 假設檢驗法標定物理量的PDF，采用Beta 分布定量刻畫沒有物理意義的、唯象參數的不確定度，評估輸入不確定度對QoIs 的影響，以增強C-J 模型的可靠性和預測精度，以期為武器型號設計、裝甲防護和土木施工管理等提供決策依據。

1 Chapman-Jouguet 理論

假設爆轟波運動方向和波陣面垂直，忽略波陣面的厚度和化學反應時間，忽略熱傳導和黏性效應以及非定常效應，則爆轟過程可簡化為圖1 所示的C-J 模型。圖1 中，ρ0、p0、U0和e0分別代表未反應物的初始密度、壓力、速度和內能，D為爆速。爆炸物以D–U0的速度進入波陣面，以D–U1的速度流出波陣面。

圖1 爆轟波在高能炸藥中的傳播示意圖Fig. 1 Sketch of propagation of detonation wave into high explosive

根據Hugoniot 守恒定律和C-J 假設，波陣面兩側物理量滿足如下關系。

(1) 質量守恒：

(2) 動量守恒：

(3) 能量方程：

(4) Chapman-Jouguet 假設：

(5) 聲速定義：

式中：v1=1/ρ1，為爆轟產物的比容；v0=1/ρ0，為未反應炸藥的比容； ρ1為爆轟產物的密度；p1為波后壓力；U1為波后介質速度；e1為爆轟產物的內能；c為聲速；S代表等熵狀態。

對于高能炸藥，p0?p1，因而假設p0=0。由于爆轟波未到達時，圖1 右側高能炸藥近似認為處于靜止狀態，進而假設U0=0。且高能炸藥反應物和產物的EOS 均使用γ 律：

由式(1)～(6)導出波后介質QoIs 的函數表達式：

式中：ρJ、pJ、UJ和cJ分別代表C-J 理論導出的爆轟產物的初始密度、壓力、速度和聲速。式(7)建立了C-J 狀態下波后介質QoIs 和物理量D、ρ0以及唯象參數γ 的函數關系式。由于波后介質物理量試驗標定困難且成本較高，因此式(7)成為預測波后介質QoIs 的實用工具。

2 爆轟不確定因素的定量刻畫

受測量技術、含能材料自身固有屬性和人類認知的局限性等因素的影響，爆轟M&S 中的物理量和唯象參數均受到噪聲污染[12-14]。ρ0的不確定度主要來源于壓制成型時炸藥顆粒的扭結、隨機結晶、空洞和裂紋的存在以及雜質的摻入；爆速D的不確定度主要來源于測量技術。即D和ρ0的不確定度與材料的固有屬性有關，即使增加樣本容量，不確定度也不會消失。

2.1 不確定物理量的表征

根據Wilkins 等[19]的理論，在爆轟計算中，合適的概率分布能合理描述含能材料密度的非均勻性。本文中，假設PBX-9502的初始密度服從對數正態分布，其中μ ?和 σ?分別為對數期望和對數標準差。與正態分布相比，對數正態分布的優勢在于符合統計規律的同時能保證物理量的非負性。首先，利用浸液法測量炸藥的密度[20]，給出PBX-9502 初始密度的統計直方圖，如圖2 所示。根據矩估計法導出試驗結果的期望μ=1.894 g/cm3和標準差σ=0.002 5 g/cm3。進而，使用Anderson-Darling 假設檢驗法，在5%顯著性水平下，ρ0接受服從對數正態分布的原假設。最后，通過曲線擬合技術，得到PBX-9502 的概率密度函數。如圖2 所示，初始密度的概率密度函數曲線與試驗數據吻合很好。

圖2 PBX-9502 初始密度的概率密度函數Fig. 2 PDF of initial density of PBX-9502

爆速是爆轟研究中唯一能使用較簡單試驗裝置和流程即可標定的物理量，這里的試驗結果來自于測時儀法[20]。假設D～，使用Anderson-Darling 假設檢驗，驗證得到：在5%顯著性水平下，D也接受服從對數正態分布的原假設。結合試驗數據，使用矩估計法，導出爆速的期望μD=7710m/s，標準差σD=321m/s。代入式(8)，計算得到=8.9503，=0.0025。使用曲線擬合技術，得到PBX-9502 爆速的概率密度函數曲線如圖3 所示。

圖3 PBX-9502 爆速的概率密度函數Fig. 3 PDF of detonation velocity of PBX-9502

2.2 不確定的唯象參數

EOS 中的參數γ 無法通過試驗直接標定。由工程經驗和專家意見，假設γ～B[2.93, 2.99, 6, 6]，其中B[a,b,m,n] 代表4 參數Beta 分布，[a,b]表示參數的取值范圍，m、n表示形狀參數。根據經驗，形狀參數取m=6，n=6。與正態分布相比，Beta 分布的優點在于取值范圍有界。與均勻分布相比，Beta 分布的PDF 曲線光滑而不間斷。γ 的PDF 的詳細信息如圖4 所示。

圖4 多方指數γ 的概率密度函數Fig. 4 PDF of polytrophic exponent γ

3 不確定度傳播方法

3.1 高維非嵌入多項式混沌

設 ξ=(ξ1,ξ2,···,ξd)T是概率空間 ( ?,F,P) 中的d維獨立標準正態隨機向量，其中Ω為樣本空間，F為定義在Ω上的σ-代數，P為概率測度。L2(Ω)為Ω上的平方可積函數空間，設u∈L2(?) 代表爆轟系統的QoI。若賦予內積：

式中： ω ∈? 代表基本事件， ρ(x)=(2π)-d/2e-xTx/2為隨機向量 ξ 的聯合概率密度函數，則L2(?) 構成Gauss-Hilbert 空間。利用Cameron-Martin 定理[17]，u可以通過 ξ 的正交多項式展開：

隨機基函數：

當d=1 時，由式(10)、(13)和Hermite-Gauss 求積公式，得到：

式中：ξr、wr分別為求積點和權重。本文中取S=6，此時，Hermite-Gauss 積分的求積點和權重如表1 所示。

表1 6 節點Hermite-Gauss 積分的求積點和權重[21]Table 1 Quadrature points and weights for Hermite-Gauss integration with six points[21]

當d＞1 時，使用全張量積Hermite-Gauss 求積公式，即：

實際應用中，式(10)通常截斷為：

3.2 Rosenblatt 變換

由第2 節知，輸入不確定度不服從標準正態分布且概率類型不一致，因此無法直接使用3.1 節所述PC 理論。本文中使用Rosenblatt 變換將一列相關隨機變量轉化成一列服從標準正態分布的獨立隨機變量。具體步驟為：

設X=(X1,X2,···,Xd)T為d維非高斯隨機向量，令FX(x1,x2,···,xd)為X的聯合累積分布函數(cumulative density function, CDF)。構造如下映射關系Y=G(X)[22]：

式中：Fxd(xd|x1,x2,···,xd-1) 為條件CDF， Φ 為標準正態變量的CDF，則Y服從獨立標準正態分布。

4 結果的UQ 分析

利用3.2 節Rosenblatt 變換技術將第2 節標定的輸入不確定度D、ρ0和γ 轉化為相互獨立的標準正態分布。然后利用3.1 節非嵌入PC 結合式(7)計算波后QoIs 的概率密度函數。具體流程如圖5 所示。

圖5 不確定分析算法流程圖Fig. 5 Flow chart of uncertainty analysis

本文中，d=3，Q=3，因此K=20。接著利用表1，構造二維全張量Hermite-Gauss 積分求積點，如圖6所示。本文中計算將會用到的前6 階Hermite 多項式的具體形式如圖7 所示。結合式(1)～(7) 和(15)～(16)，利用圖5 所示的不確定度分析流程計算出QoIs 的概率密度函數，更多信息如圖8 所示?？梢钥闯?，波后QoIs 的概率密度函數沒有明顯的對稱性，且峰度指標均大于3，處于尖峰狀態，不服從Gauss 分布。特別是波后密度，峰度指標最大，曲線最尖峭，且具有明顯的偏差，與正態分布差別很大。波后壓力的峰度系數最小，對稱性最好，形狀最接近正態分布，可用正態分布近似替代。

圖6 二維全張量積Gauss-Hermite 求積點的位置Fig. 6 Locations of Gauss-Hermite quadrature points used for two-dimensional tensor product

圖7 前6 階單變量Hermite 多項式Fig. 7 The first six orders of univariate Hermite polynomials

圖8 系統響應量的概率密度函數Fig. 8 PDF of system response quantities

此外，QoIs 的統計量，如期望和方差，可通過uα顯性表達，即：

利用式(18)計算出波后QoIs 的期望，標準差，具體內容見表2。從表2 可以看出，波后介質密度的標準差與期望的比值小于5%，而波后介質壓力的標準差與期望的比值大于5%。波后介質壓力波動較大，取值較寬，與文獻[6, 23]的“爆轟壓力測量值分散性較大”的結論相吻合。

表2 波后感興趣量的試驗和統計結果Table 2 Statistical and experimental result of QoIs at the rear of the shock wave

根據圖8 信息，可計算出95%置信水平下波后QoIs 的置信區間。同時，波后介質的壓力、密度和速度等物理量也能通過試驗標定[6,23]。數值預測結果與洛斯阿拉莫斯國家實驗室（Los Alamos National Laboratory, LANL）的試驗結果[24]作比對。比較發現，波后介質的壓力、密度和速度的試驗數據落入95%置信水平下的置信區間內，進一步確認了方法的有效性。

5 結果與討論

C-J 理論建立了波后介質壓力、速度、聲速和密度等系統響應量和初始密度、爆速以及多方指數的顯式函數關系。初始密度易于測量，爆速也容易通過較簡單的試驗標定。C-J 理論成為試驗之外的波后物理量標定手段，成本低，實用性強。

本文中，在概率框架下研究C-J 理論，結果以PDF 表征，允許QoIs 在一定范圍內變動。首先給出了輸入不確定度的定量表征，使用Anderson-Darling 假設檢驗法驗證概率假設的正確性。運用收斂速度更快的非嵌入PC 結合Rosenblatt 變換研究了輸入不確定度對波后QoIs 的影響，發現試驗結果落入QoIs 的置信區間內。這說明密度、爆速等物理量在生產、加工和測試過程中的正常擾動不影響C-J 理論的使用。同時，模型也允許多方指數在一定范圍內變動?；贑-J 理論的爆轟UQ 結果能為裝藥設計、裝甲防護等提供決策依據。方法具有普遍性，能推廣到其他狀態方程描述的爆轟UQ 研究。

下一步擬減少假設、簡化和近似的使用次數，將C-J 理論UQ 研究推廣到更加復雜的情況，如非γ 律EOS，甚至推廣到ZND 理論的UQ 研究。對于非理想爆轟，聲速面后可能存在放能。引入放能后，預測結果的差異，即不同模型導致的不確定度，也叫模型形式不確定度量化。這是國際熱點，也是不確定度量化領域的難點之一。

另外，本文中，研究沿用Oberkampf 的思路[25]，即將試驗結果看作是基準的，用以評估模擬結果的可信度。事實上，QoIs 的試驗結果也受到不確定因素的擾動，下一步擬綜合試驗不確定度和數值模擬不確定度開展研究。同時，從爆轟研究角度看，更希望從包含不確定度的某一參數實測結果中通過統計分析得到盡量多相關參數的均值和標準差等信息。若能指導實驗設計方法，用盡量少的實驗次數得到QoIs 的結果，是另一個值得探討的課題。