基于自適應和投影Wiener混沌的圓筒實驗不確定度量化*

梁 霄，王瑞利

（1.山東科技大學數學學院，山東 青島 266590；2.北京應用物理與計算數學研究所，北京 100089）

爆轟是極其復雜的物理化學過程，發生在極短的時間尺度(10～20 μs)和極小的空間尺度((1～9)×10-4m)，是爆炸力學的重點和難點。理論、實驗、數值模擬是研究爆轟的三種途徑，且這三類方法各有千秋，又相輔相成。實驗能直接反應真實物理現象，探索爆轟的發生機理。但受環境和儀器的限制，往往只能表征初始狀態和最終結果.特別是部分炸藥由于感度過高，又導致實驗人員處于危險的操作環境。數值模擬可以全時空的重復模擬爆轟的發生過程，是研究爆轟的一條安全經濟方法。同時數值模擬可以與理論、實驗相互印證，重視實驗結果，并給出了沒有實驗結果的預測[1-6]。然而由于爆轟的極端復雜性以及認知的缺陷，導致多物理爆轟數學模型特征鮮明。首先，模型結構及其運算格式高度復雜，尤其高精度數值模擬的研發和運算成本高昂。其次，目前狀態方程（EOS）和反應率方程均為唯象建模，含有大量經驗參數，沒有物理意義，根據經驗將其限定在某個范圍。最后，測量技術受外界因素的干擾，以及物質本身固有的不確定性也會導致待測物理量的隨機波動。以上原因導致真實的爆轟模型是一個眾多不確定因素耦合在一起的非線性耦合Euler方程。要發展高可信度的爆轟過程數值模擬程序，必須對炸藥起爆、描述反應區的反應率函數、狀態方程中的眾多不確定性參數進行量化。

目前，受制于爆轟模型隱藏的大量不確定度和受限于爆轟多物理過程的計算成本，多隨機因素擾動的爆轟建模與模擬的不確定度量化（uncertainty quantification，UQ）在國際上仍然是一個極具挑戰性的課題，且實質性的結果不多。然而，量化和評估不確定因素對系統輸出結果的影響，直接關系到模型的穩健性和數值模擬的可靠性。

研究爆轟模型的不確定度量化，普遍的做法是使用Monte Carlo方法(MC)。MC不需要任何假設，操作簡單，不依賴于模型的幾何結構，可以認為是系統在不確定意義下的精確解。然而，MC較慢的收斂速度(O(N-1/2))，導致采樣點數量介于103～106之間的MC方法很難應用于需要較高運行機時的高精度爆轟格式，這是MC方法在處理高維不確定爆轟模型的局限性之一。Wiener提出的多項式混沌（polynomial chaos，PC）方法在滿足特定條件下可以替代MC方法，Ghanem又將PC和有限元結合，從此PC在工程領域應用廣泛[7-9]。按照求解系數方式的不同，PC又分為嵌入多項式混沌(IPC)和非嵌入多項式混沌(NIPC)。IPC需要不斷地更改程序，而爆轟軟件的開發往往是一個團隊多年工作的結晶，因此此方法求解爆轟模型并不現實。復雜工程問題中常用的高維不確定度NIPC方法有求積法、回歸方法、采樣法等。但是鑒于高維不確定系統截斷長度和維數災難，NIPC很難直接應用于爆轟模型。

而從截斷長度和取點規模來看，上述方法會遇到如下問題。求積法遇到的一個問題是維數災難問題。簡而言之，對于一個含有20個不確定性參數的爆轟系統，每個參數選擇5個求積點，則一次積分需要運行系統520≈1014次，對于4階展開需要運行系統約1020次。對于高精度復雜多物理爆轟過程這是不現實的。至于采樣法，包含傳統的MC采樣和替代性的拉丁超立方體采樣等，目前采樣規模仍很龐大，不能直接適用于復雜高精度爆轟系統。對于回歸方法，由于爆轟計算成本較高，采樣點的個數小于切斷長度，回歸方法使用采樣將其轉化為一個欠定線性方程組，然后使用壓縮感知方法求解相應的欠定線性方程組。此方法大幅度減少了采樣的規模，Huan等[10-11]已經將其應用計算高度復雜的未來高超音速燃燒發動機數值模擬。但是這種方法至今沒有經過嚴格數學推導的誤差收斂速度和數值精度，理論上不具有完備性，使用存在未知的風險。

本文中，用自適應基函數和投影的多項式混沌方法處理含有高維不確定度的爆轟問題，并以圓筒模型為例，給出UQ結果。基本思想是：將高維隨機變量做旋轉變換，得到一個新的隨機向量組，在此隨機向量組的基礎上，重構Hermite多項式，進而利用Cameron-Martin公式重新展開，同時選取部分基函數生成相應的線性子空間。最后，將展開式在子空間上投影。 一個明顯的優勢是可以減少截斷長度，例如對于20個不確定性參數的爆轟系統，標準的4次展開PC截斷長度為(20+4)!/(20!4!)-1=10 624，通過自適應基變換和投影變換，截斷長度為(1+4)!/(1!4!)-1=4，從而在一定程度上提高了計算速度，緩解了‘維數災難’，并且此方法的誤差有精確的誤差公式。

然而，高維不確定度爆轟模型自身獨特的特征導致其還必須處理如下問題。首先，JWL狀態方程中的參數是相關的，不具備獨立同分布(independent identical distribution, IID)。其次，自適應基函數在Wiener-Hilbert空間是完備的，而基于經驗的參數未必服從正態分布。再次，以往的假設參數服從均勻分布，但是均勻分布的強間斷性，令其轉化為正態分布有一點難度。

圓筒實驗是研究爆轟波傳播驅動過程的基準實驗。在實驗中，研究人員將藥柱置于金屬圓筒內，通過多普勒類光學速度計、狹縫掃描相機獲取金屬圓筒外界面位移曲線和殼體膨脹速度，完整檢驗爆轟模型的合理性。將本文所述方法應用于圓筒實驗為基準實驗的UQ研究。也可以為后續眾多復雜物理過程的數值模擬提供條件，并能推廣到其他武器物理過程的UQ。

1 爆轟流體力學模型及不確定度

1.1 爆轟流體力學方程組

爆轟流體模型為守恒原理的Euler方程耦合化學反應的常微分方程以及復雜非線性狀態方程，在拉格朗日坐標系下的質量守恒、動量守恒、能量守恒、狀態方程和反應率方程分別為：

式中：ρ、u、E、e、p分別表示密度、速度、總能、內能與壓力，uu是指并矢張量，u·u是指內積，E=e+1/2u·u，0≤F≤1為爆轟產物的燃燒函數。

1.2 反應率方程

采用Wilkins反應率模型研究爆轟, 其形式為：

式中：F1為CJ比容燃燒函數，F2為時間燃燒函數,nb可調參數。有：

式中：v=1/ρ表示比容，v0是初始比容，vJ=γv0/(γ+1)是 CJ比容，γ是多方指數（理想氣體常數），tb是起爆時間，指爆轟波到達計算網格的時間，通過惠更斯原理計算[12]。ΔL=rbΔR/DJ，ΔR表示網格寬度，DJ是爆速，rb可調。

1.3 物態方程

爆炸產物常采用JWL狀態方程，其形式為：

式中：相對體積V=v/v0, 唯象參數A、B、R1、R2、ω可調且相關, 利用爆轟波陣面上的守恒關系以及CJ爆轟條件，得：

式中：pJ、VJ是炸藥CJ狀態下的爆壓、比容，Q是爆熱。

實際應用中，用燃燒函數F將炸藥和爆轟產物的狀態方程連接起來p=FpEOS計算反應區。

1.4 計算方法

本文的模型計算采用自主研發、具有完全知識產權的爆轟流體力學軟件LAD2D。軟件的基本算法結構如下：代碼的空間離散化基于任意多邊形非結構網格，將拉氏自相容有限體積方法及Landshoff一次黏性aL、沖擊波黏性、von Neumann-Richtmyer二次黏性aNR、子網格壓力、人工熱流等抑制非物理振蕩的方法，以及多介質滑移計算推廣到任意多邊形非結構網格，構造了一套自適應AMR任意多邊形非結構網格、流體大變形強適應的高精度有限體積格式，求解全耦合質量、動量、能量守恒方程和化學反應率方程。并采用了網格大變形處理的鄰域可變技術。時間離散化基于預設條件雙時間步。利用全顯式多步格式構造系統的隱式解，使其具有4階時間精度，并與預設條件耦合。

1.5 爆轟模型中的不確定度

關于不確定度的描述，目前并沒有統一的說法。這里不確定度有兩種類型：唯象不確定性參數和不確定性物理量。不確定度量化就是考慮兩種不確定度對系統響應量的影響。表1給出了爆轟流體力學模型的輸入不確定度。

表1 爆轟流體力學中的不確定度來源Table 1 Sources of uncertainty in detonation hydrodynamics

表中，B[α?，β?，a?，b?]代表參數?服從具有雙參數α?和β?、取值范圍在[a?，b?]的Beta分布。做線性變換Z=(?-a?)/(b?-a?)，則Z滿足標準Beta分布B[α?，β?]。由于測量的不精確性以及認知的局限性，爆轟流體力學中的唯象參數，由于沒有物理意義，因此無法通過實驗標定。這類參數的取值范圍取決于工程設計的接受程度。特別指出，以往的研究假設這類唯象參數滿足簡單易行的均勻分布，而本文中假設參數滿足Beta分布，原因在于均勻分布的密度函數的強間斷性導致其不易轉化為正態分布。此外，N(μρ,σρ2)表示初始密度ρ服從期望μρ、方差σρ2的正態分布。ρ與其余唯象參數不同，它具有明確的物理意義，產生于晶體的錯位以及炸藥凝結過程中顆粒的不均勻性。由于樣本容量大，μρ和σρ2的選取來源于統計觀測數據。炸藥選取JOB-9003, 根據實驗統計結果μρ=1.845 g/cm3,σρ=0.005 g/cm3。

圖1 不確定度的概率密度函數Fig.1 Probability density function of uncertainty

Beta分布函數的參數的確定取決于專家建議和工程經驗, 本文中參數具體取值如下：

本文中不確定度的概率密度函數如圖1所示。

2 基于自適應基函數的模型簡化

2.1 Wiener-Hilbert多項式混沌

令代表完備概率空間, 其中Ω為樣本空間，是Ω上的σ代數，P是定義在上的概率測度，本文使用 Gauss測度。θ∈Ω為基本事件。代表Ω上的分量服從標準正態分布的隨機向量，且分量獨立同分布。本文中選取d=11。令H代表Gauss空間, F(H)代表H上的σ代數。H：n：代表n次齊次∫Wiener噪聲集合。同時令L2(Ω)=L2(Ω, F (H),P)代表Ω上的∫平方可積空間，賦予內積=?f(θ)g(θ)dP(θ)，則構成一個Hilbert空間。由基礎概率論，知 =?f(θ)dP(θ)=Ef, E表示期望算子, 且滿足L2(Ω)=⊕nH：n：。

上的隨機函數，且滿足式(1)～(12)，代表密度、壓力、速度分量、總能、管壁位置等響應量。均方可積函數U(t,x,ξ)∈L2([0,T]× O ×Ω,([0,T]× O ×Ω), dt×μ×P)，其中 F ([0,T]× O ×Ω)代表[0,T]× O ×Ω 上的 σ 代數，μ是Rk上的Lebesgue測度，dt×μ×P代表([0,T]× O ×Ω)上的乘積測度，且滿足：

則根據Cameron-Martin定理[13-15],U(t,x,ξ)具有如下形式：

其中α=(α1,α2, …,αd)表示指標集, Ip={α|α1+α2+···+αd≤p}，且：

式中：Hα(ξ)代表d維 Hermite多項式。

在實際應用中，式(13)需要截斷有限次：

且根據Cameron-Martin定理，幾乎處處，。

2.2 基變換

A表示Rd上的正交變換，滿足AAT=I，I表示Rd上的單位矩陣，定義：

因此η也是H的一組基，且H：n：也可以由η生成，且令：

因此：

式中

2.3 基于投影法的模型簡化

令則U(A)(t,x,η)在VI上 的U(A,I)(t,x,η)投影定義如下：

另外，U(A)(t,x,η)在VI上的投影，又表示為：

從而：

將U限制在UI={Uγ,γ ∈ I}。

2.4 非嵌入方法求解高斯自適應基函數的系數

通過第2.3節的敘述，可以看到，模型簡化的關鍵在于A和 I 的選取。本文中使用高斯自適應方法選取A和 I。若U在高斯空間H的投影已知，則A構造如下。

首先令

式中：ei=(0, ···, 1, ···, 0)，第i個位置為 1，其余位置全為 0。

不難看出η1包含U的全部高斯統計量。η的其余分量通過Gram-Schmidt方法求解, 從而解出A。對于 I 的選取，令其包含η1的小于等于P的指標, 即 I=?1， ?1是 Ip的子集, 第i個位置不為0，其余位置全為0，因此ei是 ?1的單位向量。此時

且：

誤差為：

誤差量的期望為0，且與高斯量正交。

式(22)的系數通過非嵌入方法求解，方法如下：

式中：η(r)=(η1(r), 0, ···, 0)，η(r)和wr分別為求積點和權重，s為求積點的個數。且U(t,x,ξ)滿足式 (1)～(12),并滿足關系式 ξ=A-1η=ATη。

2.5 逆累積分布函數和Rosenblatt變換

由第2.1節知，自適應基函數理論成立的前提是{ξ1,ξ2, ···,ξd}為一列獨立同分布的標準正態隨機變量組。然而，這并不適用于本文的研究對象。由隨機變量組（見表1），{ξ1,ξ2, ···,ξn}不僅含有非正態分布隨機變量，而且即使服從正態分布也不是標準正態分布。由第1.3節式(10)～(12)知，即給定R1、R2、ω可計算對應的A、B，即A、B可由R1、R2、ω給出。因此A、B、R1、R2、ω并非相互獨立。因此，必須對隨機變量組做適當改進，方可使用第2.1～2.4節方法。

使用逆累積分布函數變換將一般Beta分布化成標準正態分布。具體步驟如下：

首先，設X、Y分別滿足累積分布函數FX(x)、FY(y)，利用概率等交換原則，令FX(x)=FY(y)，則

假設X～ N (0, 1)，則：

若Y～ N (μ,σ2), 則X=(Y-μ)/σ服從標準正態分布。

其次，使用Rosenblatt變換[16]，將一列相關隨機變量化為服從正態分布的獨立隨機變量組。利用概率等交換原則：

其中：

表示條件概率。從而：

則{Y1,Y2, ···,Yn}服從標準正態分布獨立同分布。

進而，利用Rosenblatt變換，可將不確定參數化為一列服從標準正態分布的獨立同分布隨機變量。

3 圓筒實驗的不確定度量化

圖2 圓筒實驗裝置示意圖Fig.2 Schematic diagram of cylinder test

圓筒實驗是確定炸藥爆轟產物JWL狀態方程和評估炸藥做功能力的標準化實驗，應用廣泛。實驗原理是將炸藥放入等壁厚的銅質圓筒中，從圓筒的一端將其引爆，利用高速轉鏡式掃描相機記錄筒壁在爆轟產物驅動下的膨脹過程。圓筒實驗布局見圖2，炸藥尺寸為 ? 25.0 mm×305 mm，炸藥為JOB-9003炸藥，實驗測得爆速為8 712 m/s；圓筒內徑25.0 mm，壁厚2.5 mm，材料為紫銅。狹縫掃描采用平行光后照明技術，狹縫位置距離起爆端200 mm。計算格式采用方法見第1.4節。圖3是針對JOB-9003炸藥25.0 mm圓筒實驗測試結果，給出了距起爆端200 mm狹縫處，筒壁在爆轟產物驅動下的膨脹過程，徑向壁位置與徑向壁速度隨時間的演化過程。

圖3 JOB-9003圓筒實驗結果Fig.3 Experimental results of JOB-9003 in cylinder test

利用本文所述的方法，給出了圓筒實驗的不確定度量化結果，可以從多角度觀測輸入不確定度對輸出結果特別是速度和管壁位置的影響。具體內容如圖4～8所示。圖4～5給出了管壁位置和速度隨時間變化的期望值及標準差。可以看出，爆轟波在25 μs到達管壁，并繼續向前傳播，并在30 μs附近速度到達峰值，管壁在慣性作用下繼續向前運動，同時速度峰值最大的點也是速度標準差最大的點。從圖5可以看出，本文的計算格式在強間斷處帶有明顯的波后震蕩特征，這可能也是激波位置標準差最大的部分原因。而管壁位置的期望和標準差在激波達到后呈現一定的類線性關系。

圖6給出了圓筒管壁位置和速度隨時間變化的置信區間的計算結果，區間寬度滿足3σ法則，并將置信區間涂黃。計算結果看出，爆轟波過后，置信區間寬度變大。波前標準差為0，因此區間寬度也是0。在圖7中計算結果和實驗結果對比發現，實驗樣本均落在置信區間內，符合預期。為觀測方便，圖8為圖7的局部放大，進一步驗證了觀測結論。同時，也確認了本文方法的有效性。也說明，本文中參數在特定的有限時間內可以在一定范圍內變動，均可符合精度要求。同時，進一步觀測試驗數據，發現爆轟波到達后，速度的實驗數據也有明顯的震蕩，這與標準差的震蕩吻合。說明我們的算法是可信的。

圖4 管壁位置的期望與標準差Fig.4 Expectation and standard deviation of cylindrical wall position

由于本文為含有11個不確定性參數的爆轟系統，標準的4次展開PC截斷長度為(11+4)!/(11!4!)-1=1 365，通過自適應基變換和投影變換，截斷長度為(1+4)!/(1!4!)-1=4，從而在一定程度上提高了計算速度，緩解了維數災難。

圖5 管壁速度的期望與標準差Fig.5 Expectation and standard deviation of cylindrical wall velocity

圖6 管壁位置和速度的置信區間Fig.6 Confidence intervals of position of cylindrical wall position and velocity

圖7 實驗數據和計算結果的置信區間Fig.7 Confidence intervals of experimental and simulation results

4 結論與展望

通過Wiener-Hilbert空間中基函數的旋轉變換和投影變換，明顯減少了截斷長度，從而減輕了計算量。并將此方法應用于爆轟圓筒實驗，給出了不確定度量化和傳播結果。結果以期望、標準差、置信區間表示，并與實驗數據做比較，發現實驗數據均落在置信區間內。同時，觀測到波后速度震蕩劇烈與計算結果中強間斷出數值模擬標準差達到峰值吻合。從而確認了模型的有效性。

下一步，將本文應用到其他模型并與現有的結果對比, 同時開展如下工作。

（1）本文中研究局限在Wiener-Hilbert空間，導致隨機變量必須服從正態分布。對于非正態分布變量通過等概率原則和Rosenblatt變換將其轉化成正態分布變量，這增加了計算量。而眾所周知的是，在PC理論中均勻分布對應勒讓德多項式，Beta分布對應Jacob多項式等。能否在其他完備的Banach空間中研究此自適應投影問題，從而不需要轉化，減少計算量？若把唯象參數當做認知不確定度，如何處理爆轟系統的不確定度量化？這仍是一個具有挑戰性的課題。

（2）波后震蕩與實驗數據吻合，沒必要消除，但是能否使用高精度格式改進模型算法，減少波后震蕩?需進一步提高模型的精確性。

(3)本文中主要研究參數和物理量不確定的爆轟系統的不確定度量化，然而爆轟系統的復雜性導致狀態方程和反應率方程的選擇也是不確定的，如何研究不同形式狀態方程和反應率方程對系統的影響，即模型形式不確定度的量化？這也是一個具有挑戰性的課題。

(4)僅研究炸藥狀態方程和反應率方程中的不確定參數是不夠的。在爆轟實驗中，無論是光學掃描方法還是激光干涉儀，仍然有很多不確定因素值得研究，如測量不確定度等。需要進一步深入分析實驗中的不確定度來源，量化不確定度的傳播，提高數值模擬的可信度和模型的預測能力。因此，實驗不確定度是下一步研究的課題。