脈沖X射線輻照鋁板氣化反沖仿真的不確定度量化方法

王 棟，湯文輝?，張 昆，冉憲文，晏 良

(1.國防科技大學物理系，長沙410073；2.國防科技大學體系科學系，長沙410073)

在納秒級極短脈沖寬度內，高能注量X射線輻照會導致材料比內能急劇升高，材料迎光面氣化，由此產生的熱壓分布和氣化反沖沖量(blow-off impulse,BOI)在材料內形成激波，對航天器和導彈等防護材料造成毀傷[1]。

對于脈沖X射線輻照誘導的BOI熱-力學效應，國內外從20世紀就開展了大量的理論求解和實驗研究。在理論求解方面，推導了一系列的解析計算公式，如BBAY公式[2]、Whitener公式[3]、MBBAY公式[4]及 Zhang公式[5]等。這些公式都基于一定的假設和簡化處理，主要用于快速估算。在實驗研究方面，美國國家實驗室建設了SATURN、DOUBLE EAGLE等高能注量X射線源，并開展了一系列實驗；我國在“強光一號”上開展了軟X射線輻照實驗。但上述X射線源實驗成本很高，且輻照靶材較小，對結構響應研究具有局限性。此外，我國還開展了一系列炸藥加載、柔爆索加載、電子束和離子束加載等模擬加載實驗[7-10]，但是不可避免地在實驗等效性上存在誤差。

數值仿真技術彌補了X射線源實驗條件的不足，與理論解相比，可以利用更精細的計算模型和復雜的初、邊界條件，借助超算進行大尺度仿真，發(fā)揮建模靈活性及數據采集與參數優(yōu)化等方面的優(yōu)勢。美國空軍實驗室開發(fā)的1維仿真程序PUFF-TFT[11]，考慮了次級粒子輸運引起的劑量增強效應，更精確地描述了X射線的能量沉積過程，主要用于評估X射線沉積導致的1維應變-應力響應和薄層疊加效應。美國圣地亞國家實驗室開發(fā)了CTH程序，用于建模仿真大形變和高強度沖擊下的3維復雜響應問題[12]，包括黏塑性和率相關模型，實現了自適應網格細化和大規(guī)模并行計算，有效提高了計算精度和計算效率。湯文輝等采用光滑粒子流體動力學方法對金屬鋁在脈沖X射線輻照下的BOI進行了數值模擬[13]。蔣邦海研究了某種纖維增強型復合材料的1維應變率相關動態(tài)本構模型，用于1維X射線輻照動力學的數值模擬[14]。黃霞等建立了2維正交各向異性動態(tài)彈塑性本構模型，基于有限元法開發(fā)了2維仿真程序TSHOCK2D[15]，對碳酚醛材料中的2維X射線熱擊波和BOI進行了研究。ZHANG等將各向異性本構模型推廣到3維，并在文獻[15]的基礎上開發(fā)了3維仿真程序TSHOCK3D[16]。這些仿真程序，缺少動態(tài)內存分配、動態(tài)時間步及自適應網格技術等計算優(yōu)化方法，如采用高精度網格，仿真計算將耗時較大。

在對脈沖X射線輻照動力學過程的建模與模擬中，由于存在模型及材料物性參數測量不確定度等情況，仿真輸入參數中存在隨機不確定性和認知不確定性，通過計算模型的傳導，導致材料響應和結構響應的估計結果的不確定性，最終會導致氣化反沖仿真結果的不確定性，對數值仿真結果的可靠性或可信度評估帶來挑戰(zhàn)。如何量化數值仿真中由于輸入參數不確定性導致的響應量的不確定性是一個值得探討的問題。

針對MC方法收斂性低的問題，一種解決方法是使用擬(MC)方法，該方法的不確定度近似為O(1/N)量級，但仍有一定的局限性；另一種解決方法是用代理模型(也稱元模型)，近似原仿真模型，可顯著減少計算時間?，F有的代理模型眾多，如響應面法、徑向基函數法及樣條插值法等[19]。與其他方法相比，Kriging代理模型具有回歸和插值的雙重特性及良好的非線性估計性能，能給出預測點處的估計方差，對開展不確定量化和模型精度評估具有一定優(yōu)勢，應用較廣。

但是，本文仿真的高精度樣本較為稀缺，計算耗時僅能承受個位次數的仿真，不足以構造足夠精確的Kriging模型來解決12維的問題。因此，本文通過網格收斂性分析，選取一組具有均衡時效比的低精度網格來獲取低精度樣本，將高、低精度樣本通過訓練Co-Kriging模型進行數據融合，有效解決了數值仿真樣本點獲取成本較高帶來的代理模型精度較低，甚至無法建立代理模型的問題；同時，提出了一種基于留一法交叉驗證的序貫試驗設計方法，提高了采樣效率，確保了代理模型精度；最后，基于TSHOCK3D實現了X射線輻照鋁板氣化反沖仿真，并進行了不確定度量化分析。

1BOI的有限元仿真

X射線的輻射功率通常采用Planck黑體輻射模型：

(1)

其中，c1，c2分別為第一、第二輻射常量；λ為波長；T為黑體溫度。在數值上，將λ離散，并截斷上下限。歸一化處理后，對于波長處于Δλi內的光子，能量占入射能量的份額ωi(λ)為

(2)

X射線能量沉積計算中主要考慮光電效應和Compton散射效應。查表可得Δλi內鋁對光子的質量吸收系數μi(λ)。假定X射線脈沖波形是矩形波，在1維條件下，X射線能注量Φ(x)在材料內隨入射深度x的變化關系為

(3)

其中，n為離散個數；ρ0為材料的初始密度；Φ0為初始入射能通量。對于3維條件下的能量沉積計算，采用矢量面積分方法進行簡化計算[20]。

BOI的數值計算方法有2種。一種是通過計算迎光面處的壓力對時間積分獲得；另一種是通過計算氣化單元的動量之和獲得。文獻研究認為，氣化單元尺寸的急劇膨脹，使得單元兩端節(jié)點處的速度和密度差較大，進而導致數值計算結果的不確定度增加[21]。因此采用第一種方法計算更為準確。

更新拉格朗日方法的控制方程組為

(5)

采用彈塑性流體力學本構，塑性屈服采用von-Mises屈服準則，斷裂采用最大拉應力準則，物態(tài)方程采用修正的各向異性PUFF物態(tài)方程[22]。

2不確定度量化方法

近年來，不確定性量化理論快速發(fā)展，但目前還沒有統一、標準的數學框架。基于樣本的不確定性量化分析有4個基本步驟[23]：1)通過抽象物理模型，構造數學模型和計算模型，建立輸入變量x和目標響應量y的數值映射關系；2)根據專家經驗或統計數據，利用概率密度分布函數描述輸入變量x的認知不確定度或隨機不確定度；3)從輸入變量x到目標響應量y進行不確定度傳遞；4)給出不確定性分析結果，如目標響應量的經驗概率密度分布。此分析過程一般利用MC方法，通過對隨機輸入參數進行概率抽樣和確定性仿真，完成不確度傳遞，并對目標響應量y進行概率統計分析。

本文基于樣本的不確定性量化分析框架，對不確定度量化過程進行細化和改進，利用Co-Kriging代理模型來融合兩類仿真精度數據點，進行MC仿真，不確定度量化流程如圖1所示。

圖1不確定度量化流程圖Fig.1Flow chart of uncertainty quantification

2.1Kriging模型的構造

常見的Kriging代理模型有Simple Kriging (SK) 模型，Ordinary Kriging (OK) 模型，Universal Kriging (UK) 模型和Blind Kriging (BK) 模型，主要差異在于假設的回歸函數的形式不同。選擇回歸函數是Kriging建模的難點，目前主要通過經驗或試算來選擇。工程問題一般假設回歸函數為常數，采用OK模型。對于UK模型，基函數選擇不當會降低代理模型的估計精度。BK模型通過一種貝葉斯特征選擇技術，從基函數備選集中選擇有顯著影響的基函數加入回歸函數中，需要利用更多的樣本點進行訓練，但其在某些測試函數中的性能與OK模型接近[24]。因此，對于計算參數維度高及非線性程度高的工程問題，綜合考慮計算時間和擬合精度，OK模型是較優(yōu)的選擇[25]。但是，對于樣本點抽樣困難的情況，OK模型的預測精度較差。O’Hagan等將Co-Kriging (CK) 模型引入工程仿真計算領域，利用易抽樣的低可信度樣本來輔助抽樣困難的高可信度樣本構建代理模型，提高了模型整體預測精度[26]。

綜上，由于本問題的高精度仿真計算時間較長，本文構造具有較高保真度的代理模型予以解決。首先，采用拉丁超立方設計，計算得到一定數量的高精度仿真樣本；其次，基于OK模型進行序貫試驗設計，計算得到低精度仿真樣本；最后，基于CK模型融合不同精度的仿真樣本，實現代理模型構造效率和全局預測精度的平衡。

給定樣本集X={x1,x2,…，xn}，x∈d對應的響應量集記為y={y1,y2,…，yn}，y∈,輸入輸出樣本的維數分別是d維和1維，樣本數為n。為了消除量綱的影響，提高Kriging模型的精度和魯棒性，對樣本集和響應量集進行歸一化預處理。

通常，Kriging模型包含插值和回歸兩重屬性，其構造包括2步，首先基于樣本集構造回歸函數(趨勢函數)，然后基于殘差構造高斯過程Z(x)。OK模型假設回歸函數為常數，可以寫為如下形式：

Y(x)=μ+Z(x)

(6)

其中，μ為常值回歸函數；Z(x)為高斯過程函數，滿足下面統計學特征：

其中，R(Θ,x,x′)為任意兩個樣本點x和x′之間的相關函數模型，包含一組超參數Θ={θ1,θ2，…,θn}。對工程類問題，Matérn相關函數具有更好的擬合效果[27]，這里采用Matérn32相關函數，其形式為

(8)

(10)

其中，r(x)是n維向量,第i個元素是待估樣本點x和已知樣本點x(i)的相關函數R(Θ,x,x(i))。

(12)

(13)

(14)

其中，(σc,rc,Rc)和(σd,rd,Rd)分別由OK模型Yc(x)和Yd(x)進行極大似然估計得到；回歸系數β從Yd(x)的極大似然估計中獲得。

2.2序貫試驗設計

Kriging代理模型的預測精度很大程度上取決于訓練集的數量。通過試驗設計可以用盡量少的試驗點得到盡可能多的模型信息，從而高效地構建Kriging模型。Kriging模型通常配合拉丁超立方抽樣(Latin hypercube sampling, LHS)，用一次性采樣點集訓練Kriging模型。但是一次性采樣會造成樣本點的冗余或不足，導致計算資源浪費或達不到要求的代理模型精度。

序貫設計(sequential design, SD)可以有效解決上述問題。其基本思想是用一組數量較少的初始訓練樣本點為基礎，根據已有樣本點信息和已構造的近似模型信息對精確模型的定義域進行序貫采樣，有針對性地補充精確模型信息，以逐步提高近似模型的全局或局部預測精度，直到滿足預設要求的精度[28]。序貫試驗設計的核心是根據試驗設計目標，選擇合適的序貫加點準則，本質上是利用某個目標函數來衡量新樣本點對模型精度的改善，將樣本點的選取問題轉化為最優(yōu)化問題。重復上述過程直到滿足終止條件停止采樣。

本文使用應用較廣的“最大方差”加點準則[29]，最大程度減小Kriging模型的全局預測方差，選擇參數空間內使得預測方差s2(x)最大的點xnew作為下一個采樣點：

(15)

(16)

在實際計算中，隨著樣本點增多，LOOCV方法帶來的開銷逐步增大。文獻[31]提出了一種近似計算公式，在留一交叉驗證時保持訓練樣本集不變，大大提高了計算效率。近似計算公式為

(17)

其中，Hi,:,H:,i分別代表矩陣的第i行和第i列元素，Hii代表對角線上元素。

最后，當σCVPE的連續(xù)相對改善(successive relative improvement, SRI)σSRI低于某個閾值時，采樣終止。σSRI的定義為

(18)

綜上，本文提出的序貫試驗設計算法如下：

Step 1：對于d維輸入參數，基于最大化最小距離準則進行初始LHS抽樣n0=2d+1個點，并對樣本點進行仿真，得到初始樣本集X(n0)和響應量集Y(n0)。

Step 2：訓練OK模型，計算初始σCVPE(n0)。

Step 3：選擇xnew作為下次仿真計算樣本點xn0+1，更新OK模型，并計算σCVPE(n0+1)。

Step 4：計算連續(xù)相對改善σSRI(n0)。若樣本數n≤nmin，回到Step3。

Step 5：若σSRI(n0)連續(xù)4次低于閾值σ0，則停止。否則回到Step3。

需要注意的是，為防止前期σSRI的偶然性波動，需設定一個最低迭代次數nmin=4d。

序貫試驗設計流程如圖2所示。

圖2序貫試驗設計流程圖Fig.2Flow chart of sequential experiment design

3算例與分析

利用有限元程序TSHOCK3D進行算例仿真和結果分析，使用Ansys前處理生成正交六面體單元網格，并且假設X射線垂直靶面入射，有限元網格劃分以及仿真場景，如圖3所示。

圖3有限元網格劃分以及仿真場景示意圖Fig.3Schematic diagram of finite element griddivision and simulation scene

3.1輸入參數不確定度分布

關于輸入參數的不確定度分布，本文僅根據專家經驗給出粗略估計。對于鋁的10個物性參數，采用3σ準則給出標準偏差的估計，出于簡化問題考慮，取3σi=0.1μi，即假設各參數以99.73%的概率落在μi±3σi區(qū)間內。對與輻射源有關的X射線復合譜及能注量，則基于實際X射線復合譜及能注量的等效簡化假設選取值，并根據專家經驗假設正態(tài)分布，根據3σi=0.1μi給出標準偏差的估計。

針對這一組參數x0開展氣化反沖仿真的不確定度量化，具體各參數的取值及不確定度分布按上述估計給出，如表1所列。

表1各參數的取值及不確定度分布Tab.1The distribution of random parameters

3.2網格收斂性分析

進行網格收斂性分析，一是選擇1組高精度網格，近似地滿足計算網格無關性條件；二是選擇1組合適的低精度網格來開展多精度數據的融合，在代理模型精度和樣本計算開銷上尋求平衡。算例在0.1 m×0.2 m×0.2 m的靶板區(qū)域內進行仿真，考慮到對稱性邊界只計算1/4模型。為避免邊側稀疏波的干擾，y，z方向長度為x方向長度的4倍。取靶板迎光面和x軸交匯點(原點O)為應力監(jiān)測點，對時間積分得到BOI。采用中間加密法細化網格，對得到的各組網格進行仿真，給出鋁板BOI隨網格數量變化關系，即鋁板BOI的收斂性曲線，如圖4所示。

圖4鋁板BOI的網格收斂性曲線Fig.4Grid convergence curve of blow-off impulse

由圖4可見，隨著網格密度增加，CPU計算時長快速增加，BOI逐漸收斂。原因是軟X射線的能量沉積曲線在迎光面附近急劇下降，x方向網格精度極大地影響了能量沉積離散的計算精度。隨著網格精度的提高，離散單元的能量沉積逐漸收斂，最終與連續(xù)能量沉積曲線對應。由于程序采用均勻網格，提高網格精度將帶來CPU計算時長的增加，降低邊際效益。因此，綜合考慮時效比，對照網格收斂性曲線，選定BOI收斂速度趨緩處對應網格作為低精度網格，記作Lc。將BOI基本達到收斂處對應網格作為高精度網格，記作Le。取產生BOI的時刻為0時刻，對比兩類網格仿真的BOI隨時間的變化關系，如圖5所示。

圖5采用Le和Lc精度網格的鋁板BOI隨時間的變化關系Fig.5Comparation of Le and Lc grid simulationof blow-off impulse induced by pulsedX-ray irraditation on aluminum

3.3序貫試驗設計和CK建模

在高精度樣本集的選取上，基于最大化最小距離準則對隨機參數向量x0進行一次性的LHS采樣。由于計算耗時的限制，設計6個樣本點進行Le精度仿真。在低精度樣本集的選取上，采用基于最大化最小距離準則的LHS設計初始25個樣本點，并進行Lc精度仿真。選擇t=0.5 μs時刻(平臺期)的BOI作為響應量，訓練初始OK模型，按照序貫采樣算法繼續(xù)采樣92次，直至達到停止條件。此時，迭代訓練后的OK模型σCVPE=0.003 4。σSRI的收斂趨勢如圖6所示。

圖6Lc精度網格仿真σSRI收斂趨勢圖Fig.6Convergence trend of Lc grid simulation

用117個Lc精度仿真結果和6個Le精度仿真結果中t=0.5 μs時刻的BOI作為響應量，訓練得到CK-1模型，相應σCVPE=0.009 3。同時，選擇t=0.05 μs時刻(上升期)的BOI作為響應量，訓練得到t=0.05 μs時刻的CK-2模型，相應σCVPE=0.009 0，可以看到CK-2具有與CK-1接近的模型預測精度。

3.4響應量BOI的不確定度量化

核密度估計 (kernel density estimation, KDE) 是在概率論中用來估計未知參數的密度函數，屬于非參數檢驗方法之一。基于已經獲取的MC仿真樣本，本文采用Epanechnikov核函數對響應量總體的概率密度進行估計，得到近似的概率密度函數，即經驗概率密度函數 (empirical probability density function, EPDF)。利用上述兩個時刻數據訓練的CK-1和CK-2模型，分別進行MC仿真1×106次。分別對兩個時刻的樣本點集進行核密度估計，得到t=0.5 μs和t=0.05 μs時刻的BOI的期望和標準差的經驗概率密度函數及其95%置信區(qū)間，如圖7所示。

(a)t=0.05 μs

(b)t=0.05 μs

(c)t=0.5 μs

(d)t=0.5 μs

由圖7可見，在t=0.05 μs時刻，根據MC仿真得到BOI的95%置信區(qū)間為111.40～117.29 g·cm·μs-1，實際Le精度仿真結果Ie=113.07 g·cm·μs-1落在95%置信區(qū)間內。在t=0.5 μs時刻，根據MC仿真得到BOI的95%置信區(qū)間為263.29～280.94 g·cm·μs-1，實際Le精度仿真結果Ie=265.04 g·cm·μs-1落在95%置信區(qū)間內。因此，認為以高精度網格Le采用輸入參數組x0仿真得到的BOI在5%顯著性水平下是可信的。

此外，基于建立的CK模型，還可實現±10%參數空間內的BOI的快速預估。下面重新設計6組樣本點，對比Le精度仿真結果和CK模型預測結果。其中，前6組用CK-1模型進行估計，后6組用CK-2模型進行估計。對比TSHOCK3D仿真結果，預測的平均相對偏差在5%以下，如表2所示。

表2CK模型對BOI的快速評估及相對偏差Tab.2CK model fast prediction for BOI and its relative error

4結論

1)針對脈沖X射線輻照鋁靶BOI數值仿真輸入參數的不確定性問題，本文提出了一種改進的MC不確定度量化方法。通過加入序貫試驗設計采集樣本，及采用CK模型進行數據融合，有效提高了不確定度量化的效率。給出了1 keV軟X射線脈沖輻照下鋁板BOI仿真結果的經驗概率密度分布及其95%置信區(qū)間。計算結果表明，基于高精度網格的仿真結果在5%顯著性水平下是可信的。

2)通過新的高精度仿真樣本驗證，本文建立的Kriging代理模型，平均估計相對偏差在5%以下，可用于x0(1±10%)參數空間內的BOI的快速預估。

3)下一步將根據相關實驗或數值計算及其他不確定度信息，更準確地給出分布參數的估計。同時，基于建立的Kriging模型進行輸入參數的敏感性分析，如采用Sobol法給出對仿真結果影響較大的輸入參數。