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

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

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

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

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

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

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

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

1BOI的有限元仿真

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

(1)

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

(2)

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

(3)

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

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

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

(5)

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

2不確定度量化方法

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

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

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

2.1Kriging模型的構(gòu)造

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

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

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

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

Y(x)=μ+Z(x)

(6)

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

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

(8)

(10)

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

(12)

(13)

(14)

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

2.2序貫試驗設(shè)計

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

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

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

(15)

(16)

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

(17)

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

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

(18)

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

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

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

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

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

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

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

序貫試驗設(shè)計流程如圖2所示。

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

3算例與分析

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

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

3.1輸入?yún)?shù)不確定度分布

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

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

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

3.2網(wǎng)格收斂性分析

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

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

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

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

3.3序貫試驗設(shè)計和CK建模

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

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

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

3.4響應(yīng)量BOI的不確定度量化

核密度估計 (kernel density estimation, KDE) 是在概率論中用來估計未知參數(shù)的密度函數(shù)，屬于非參數(shù)檢驗方法之一?；谝呀?jīng)獲取的MC仿真樣本，本文采用Epanechnikov核函數(shù)對響應(yīng)量總體的概率密度進行估計，得到近似的概率密度函數(shù)，即經(jīng)驗概率密度函數(shù) (empirical probability density function, EPDF)。利用上述兩個時刻數(shù)據(jù)訓(xùn)練的CK-1和CK-2模型，分別進行MC仿真1×106次。分別對兩個時刻的樣本點集進行核密度估計，得到t=0.5 μs和t=0.05 μs時刻的BOI的期望和標準差的經(jīng)驗概率密度函數(shù)及其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時刻，根據(jù)MC仿真得到BOI的95%置信區(qū)間為111.40～117.29 g·cm·μs-1，實際Le精度仿真結(jié)果Ie=113.07 g·cm·μs-1落在95%置信區(qū)間內(nèi)。在t=0.5 μs時刻，根據(jù)MC仿真得到BOI的95%置信區(qū)間為263.29～280.94 g·cm·μs-1，實際Le精度仿真結(jié)果Ie=265.04 g·cm·μs-1落在95%置信區(qū)間內(nèi)。因此，認為以高精度網(wǎng)格Le采用輸入?yún)?shù)組x0仿真得到的BOI在5%顯著性水平下是可信的。

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

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

4結(jié)論

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

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

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