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

梁 霄，王瑞利

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

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

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

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

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

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

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

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

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

1.1 爆轟流體力學方程組

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

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

1.2 反應(yīng)率方程

采用Wilkins反應(yīng)率模型研究爆轟, 其形式為：

式中：F1為CJ比容燃燒函數(shù)，F(xiàn)2為時間燃燒函數(shù),nb可調(diào)參數(shù)。有：

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

1.3 物態(tài)方程

爆炸產(chǎn)物常采用JWL狀態(tài)方程，其形式為：

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

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

實際應(yīng)用中，用燃燒函數(shù)F將炸藥和爆轟產(chǎn)物的狀態(tài)方程連接起來p=FpEOS計算反應(yīng)區(qū)。

1.4 計算方法

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

1.5 爆轟模型中的不確定度

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

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

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

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

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

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

2 基于自適應(yīng)基函數(shù)的模型簡化

2.1 Wiener-Hilbert多項式混沌

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

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

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

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

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

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

且根據(jù)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 非嵌入方法求解高斯自適應(yīng)基函數(shù)的系數(shù)

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

首先令

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

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

且：

誤差為：

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

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

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

2.5 逆累積分布函數(shù)和Rosenblatt變換

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

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

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

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

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

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

其中：

表示條件概率。從而：

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

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

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

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

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

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

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

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

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

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

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

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

圖7 實驗數(shù)據(jù)和計算結(jié)果的置信區(qū)間Fig.7 Confidence intervals of experimental and simulation results

4 結(jié)論與展望

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

下一步，將本文應(yīng)用到其他模型并與現(xiàn)有的結(jié)果對比, 同時開展如下工作。

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

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

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

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