梁 霄，王瑞利

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

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

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

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

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

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

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

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

1 爆轟流體力學(xué)模型及不確定度

1.1 爆轟流體力學(xué)方程組

爆轟流體模型為守恒原理的Euler方程耦合化學(xué)反應(yīng)的常微分方程以及復(fù)雜非線性狀態(tài)方程，在拉格朗日坐標(biāo)系下的質(zhì)量守恒、動(dòng)量守恒、能量守恒、狀態(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í)間燃燒函數(shù),nb可調(diào)參數(shù)。有：

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

1.3 物態(tài)方程

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

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

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

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

1.4 計(jì)算方法

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

1.5 爆轟模型中的不確定度

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

表1 爆轟流體力學(xué)中的不確定度來(lái)源Table 1 Sources of uncertainty in detonation hydrodynamics

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

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

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

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

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

2.1 Wiener-Hilbert多項(xiàng)式混沌

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

上的隨機(jī)函數(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測(cè)度，dt×μ×P代表([0,T]× O ×Ω)上的乘積測(cè)度，且滿足：

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

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

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

在實(shí)際應(yīng)用中，式(13)需要截?cái)嘤邢薮危?/p>

且根據(jù)Cameron-Martin定理，幾乎處處，。

2.2 基變換

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

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

因此：

式中

2.3 基于投影法的模型簡(jiǎn)化

令則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ù)

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

首先令

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

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

且：

誤差為：

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

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

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

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

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

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

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

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

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

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

其中：

表示條件概率。從而：

則{Y1,Y2, ···,Yn}服從標(biāo)準(zhǔn)正態(tài)分布獨(dú)立同分布。

進(jìn)而，利用Rosenblatt變換，可將不確定參數(shù)化為一列服從標(biāo)準(zhǔn)正態(tài)分布的獨(dú)立同分布隨機(jī)變量。

3 圓筒實(shí)驗(yàn)的不確定度量化

圖2 圓筒實(shí)驗(yàn)裝置示意圖Fig.2 Schematic diagram of cylinder test

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

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

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

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

圖4 管壁位置的期望與標(biāo)準(zhǔn)差Fig.4 Expectation and standard deviation of cylindrical wall position

由于本文為含有11個(gè)不確定性參數(shù)的爆轟系統(tǒng)，標(biāo)準(zhǔn)的4次展開PC截?cái)嚅L(zhǎng)度為(11+4)!/(11!4!)-1=1 365，通過(guò)自適應(yīng)基變換和投影變換，截?cái)嚅L(zhǎng)度為(1+4)!/(1!4!)-1=4，從而在一定程度上提高了計(jì)算速度，緩解了維數(shù)災(zāi)難。

圖5 管壁速度的期望與標(biāo)準(zhǔn)差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í)驗(yàn)數(shù)據(jù)和計(jì)算結(jié)果的置信區(qū)間Fig.7 Confidence intervals of experimental and simulation results

4 結(jié)論與展望

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

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

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

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

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

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