克里金模型用于預(yù)測(cè)聚能射流速度的方法研究

雷婧宇，王利俠，寇鵬飛，薛 標(biāo)

(西安近代化學(xué)研究所，陜西 西安 710065)

引 言

聚能裝藥是目前及未來(lái)用以打擊裝甲目標(biāo)最有效的戰(zhàn)斗部形式之一[1]，它依靠炸藥對(duì)藥型罩進(jìn)行爆炸加載，從而形成定向、高速射流，以此穿透并擊毀目標(biāo)。要使聚能裝藥對(duì)目標(biāo)具有較高的打擊能力，所形成的射流必須具備速度高、連續(xù)性好等特點(diǎn)。由于聚能裝藥爆炸形成射流的過(guò)程非常復(fù)雜，而獲得射流速度需要求解非線性函數(shù)，因而往往難以獲得其精確的解析值。另外，作為有效近似手段的數(shù)值模擬法也會(huì)受到材料本構(gòu)模型及所用動(dòng)態(tài)參數(shù)等因素的制約[2]。

聚能射流的速度一般可通過(guò)定常理想不可壓縮流體力學(xué)理論模型來(lái)估計(jì)[3]，但預(yù)測(cè)結(jié)果往往過(guò)高，且不能反映射流速度梯度及射流的伸長(zhǎng)，而后發(fā)展的PER準(zhǔn)定常理論[4]，可以進(jìn)一步預(yù)測(cè)射流的速度梯度和長(zhǎng)度。

一些學(xué)者建立了理論模型或經(jīng)驗(yàn)?zāi)Ｐ?，以快速預(yù)測(cè)聚能裝藥的部分響應(yīng)量[5-8]。黃正祥等[9]應(yīng)用灰色預(yù)測(cè)模型，對(duì)聚能射流侵徹陶瓷裝甲的侵徹深度進(jìn)行了預(yù)估。近年來(lái)，神經(jīng)網(wǎng)絡(luò)逐漸被應(yīng)用到聚能裝藥的射流速度預(yù)測(cè)中[10-11]，其預(yù)測(cè)結(jié)果較為準(zhǔn)確，且計(jì)算速度相比數(shù)值模擬更快，能夠?qū)崿F(xiàn)對(duì)聚能射流速度的快速計(jì)算模擬。相比理論模型與灰色預(yù)測(cè)模型，神經(jīng)網(wǎng)絡(luò)所建立的模型靈活度以及精準(zhǔn)度更高，且模型訓(xùn)練的過(guò)程不需要復(fù)雜的理論推導(dǎo)過(guò)程，具有較好的發(fā)展前景。然而，神經(jīng)網(wǎng)絡(luò)所建立的模型往往比較龐大，且其神經(jīng)元的個(gè)數(shù)以及網(wǎng)絡(luò)層數(shù)往往難以確定。

克里金代理模型由南非采礦工程師Krige提出[12]，后來(lái)發(fā)展為統(tǒng)計(jì)學(xué)的一類插值方法?？死锝鹉Ｐ屯ㄟ^(guò)訓(xùn)練樣本與其函數(shù)值之間的關(guān)系，來(lái)預(yù)測(cè)待測(cè)樣本的函數(shù)值，其優(yōu)勢(shì)在于不僅能夠給出預(yù)測(cè)值，而且能夠提供預(yù)測(cè)方差，以評(píng)估模型預(yù)測(cè)的不確定性?？死锝鹉Ｐ筒粌H靈活性高，且能夠很好地平衡其計(jì)算精度與計(jì)算效率，因而被廣泛地用以代理隱式函數(shù)的模型[13-15]。

為了快速預(yù)測(cè)聚能裝藥的射流速度，本研究以Matlab為平臺(tái)構(gòu)建了自適應(yīng)克里金代理模型，以炸藥爆速、藥型罩錐角和藥型罩壁厚作為輸入，以聚能射流速度作為輸出，通過(guò)自學(xué)習(xí)方法，獲得了預(yù)測(cè)精度較高的克里金模型。該模型可以實(shí)現(xiàn)對(duì)金屬射流速度快速、準(zhǔn)確地預(yù)測(cè)，以用于指導(dǎo)聚能裝藥的理論設(shè)計(jì)。應(yīng)用LS-DYNA有限元軟件[16-19]對(duì)所設(shè)計(jì)的聚能裝藥結(jié)構(gòu)進(jìn)行模擬，得到其射流速度，將相應(yīng)參數(shù)下獲得的數(shù)值計(jì)算結(jié)果與克里金模型的預(yù)測(cè)結(jié)果進(jìn)行誤差分析，以驗(yàn)證克里金代理模型的預(yù)測(cè)精度。

1 克里金代理模型的基本理論

(1)

其中h(x)=[h1(x),h2(x),…,hp(x)]T表示關(guān)于x的回歸模型，可取0階、1階、2階多項(xiàng)式；β=[β1,β2,…,βp]T為回歸參數(shù)向量；p的取值如下：

(2)

z(x)是均值為0的穩(wěn)態(tài)高斯過(guò)程，其協(xié)方差為：

cov(x(k1),x(k2))=σ2Rθ(x(k1),x(k2))

(3)

式中：x(k1)和x(k2)表示樣本的兩個(gè)實(shí)現(xiàn)值；σ2和Rθ(x(k1),x(k2))分別為穩(wěn)態(tài)高斯過(guò)程z(x)的方差和相關(guān)函數(shù)。本研究使用如下的相關(guān)函數(shù)：

(4)

其中θi(i=1,2,…,n)是參數(shù)。

根據(jù)克里金模型的理論，可得回歸參數(shù)向量β和穩(wěn)態(tài)高斯過(guò)程z(x)的方差σ2如下：

(5)

(6)

式中：F為回歸矩陣，其元素為Fji=hi(x*(j))(j=1,2,…,N*;i=1,2,…,n)；R為相關(guān)矩陣，其元素為Rl1l2=Rθ(x*(l1),x*(l2))(l1,l2=1,2,…,N*)。

(7)

(8)

其中r(x)=[Rθ(x,x*(1)),Rθ(x,x*(2)),…,Rθ(x,x*(N*))]T且u(x)=FTR-1r(x)-h(x)。

2 聚能裝藥射流速度的數(shù)值計(jì)算

本研究利用LS-DYNA有限元軟件，對(duì)錐角2α分別為46°、50°、60°，罩壁厚δ為1.1、1.2、1.3mm，使用B炸藥、RDX基炸藥和HMX基炸藥的27個(gè)聚能裝藥模型形成的射流進(jìn)行模擬計(jì)算，獲得射流形態(tài)表征參數(shù)以及射流速度。

2.1 聚能裝藥結(jié)構(gòu)

采用直徑Ф50mm的聚能裝藥，其結(jié)構(gòu)由藥型罩、炸藥裝藥、傳爆起爆裝置三部分組成，藥型罩采用紫銅材料，裝藥結(jié)構(gòu)長(zhǎng)徑比大于1.5，如圖1所示。

圖1 聚能裝藥結(jié)構(gòu)Fig.1 Schematic diagram of the shaped charge

2.2 計(jì)算模型建立及參數(shù)

利用LS-DYNA有限元軟件，針對(duì)射流成型問(wèn)題的特性，采用Euler網(wǎng)格建立包含藥型罩、炸藥和空氣介質(zhì)三部分的數(shù)值模型，采用多物質(zhì)ALE算法來(lái)模擬炸藥爆轟、藥型罩壓垮及射流形成過(guò)程。其中，炸藥裝藥均采用關(guān)鍵字*MAT_HIGH_EXPLOSIVE_BURN材料模型和JWL狀態(tài)方程進(jìn)行描述，B炸藥、HMX基炸藥以及RDX基炸藥裝藥的密度分別為1.72、1.81和1.70g/cm3，其爆轟性能參數(shù)和JWL狀態(tài)方程參數(shù)見(jiàn)表1。紫銅藥型罩采用彈塑性材料模型和Gruneisen狀態(tài)方程描述[20]?？諝饨橘|(zhì)用NULL材料模型以及LINEAR_POLYNOMIAL狀態(tài)方程，空氣密度取1.293×10-3g/cm3。

表1 炸藥爆炸性能參數(shù)Table 1 Parameters of the explosives

2.3 數(shù)值計(jì)算結(jié)果

以1.1mm厚度藥型罩為示例，分別采用B炸藥、HMX基炸藥、RDX基炸藥為裝藥，不同錐角藥型罩在起爆后40μs的模擬結(jié)果分別如圖2所示，可以看到，對(duì)于同一種炸藥裝藥，起爆后的射流伸長(zhǎng)隨錐角的增大而減小；炸藥能量提高后，射流長(zhǎng)度明顯增加。表2為數(shù)值計(jì)算得到的射流性能參數(shù)。

圖2 炸藥在不同罩錐角下形成的射流計(jì)算形態(tài)Fig.2 Forms of jet produced by the explosives under different cone angles

表2 炸藥不同錐角時(shí)射流計(jì)算參數(shù)Table 2 Jet parameters of the explosives under different cone angles

由上述計(jì)算結(jié)果可以看出，在同一種炸藥裝藥中，錐角不變的情況下，隨著藥型罩厚度增大，射流頭部速度降低，頭部直徑、尾部直徑增大，有效長(zhǎng)度減??；藥型罩厚度不變的情況下，隨著錐角的增大，射流頭部速度降低，頭部直徑和尾部直徑增大，有效長(zhǎng)度減小。對(duì)比3種炸藥裝藥的結(jié)果，炸藥能量的提高可以增大藥型罩的爆炸載荷，從而提高射流速度和長(zhǎng)度。

2.4 自適應(yīng)克里金模型預(yù)測(cè)射流速度

表征射流性能的參數(shù)包括射流速度、長(zhǎng)度以及速度梯度等都可以通過(guò)數(shù)值模擬計(jì)算來(lái)獲取，但其計(jì)算時(shí)間成本仍然較高。為了快速預(yù)測(cè)聚能裝藥在不同炸藥和不同藥型罩結(jié)構(gòu)下的性能參數(shù)，以射流頭部速度為例，以炸藥爆速、藥型罩錐角、藥型罩壁厚作為輸入，射流頭部速度為輸出，建立克里金模型，實(shí)現(xiàn)射流頭部速度的準(zhǔn)確預(yù)測(cè)。此外，為了實(shí)現(xiàn)樣本的有效利用，采用克里金標(biāo)準(zhǔn)差作為學(xué)習(xí)函數(shù)，在每次迭代時(shí)自動(dòng)加入最優(yōu)訓(xùn)練樣本，以自適應(yīng)更新克里金模型。該模型預(yù)測(cè)射流頭部速度的流程如圖3所示，具體過(guò)程如下：

(1)以炸藥爆速、藥型罩錐角和藥型罩壁厚作為輸入，記為xi(i=1,2,…,N)，N為樣本個(gè)數(shù)(本研究中N=27)，建立樣本池S(x);

(3)判斷克里金模型的收斂性。計(jì)算樣本池中所有樣本的預(yù)測(cè)值μgk(xi)以及克里金標(biāo)準(zhǔn)差σgk(xi)，找到標(biāo)準(zhǔn)差的最大值，如公式(9)所示：

(9)

如果σmax<σ*(本研究中取σ*=10-3)，則克里金模型收斂，μgk(xi)即為樣本xi的克里金預(yù)測(cè)值；

(4)如果σmax≥σ*，利用學(xué)習(xí)函數(shù)選取下一個(gè)樣本，加入訓(xùn)練集T，如公式(10)所示：

(10)

圖3 自適應(yīng)克里金模型的流程圖Fig.3 The flow chart of the adaptive Kriging model

本研究在27組樣本中隨機(jī)選取5組作為初始訓(xùn)練樣本，并利用LS-DYNA計(jì)算其射流頭部速度。在訓(xùn)練克里金模型的過(guò)程中，根據(jù)公式(10)選取新的訓(xùn)練樣本，當(dāng)27組樣本中的最大克里金標(biāo)準(zhǔn)差為5.38×10-4時(shí)，克里金模型收斂，則停止計(jì)算，最終訓(xùn)練樣本的個(gè)數(shù)為11。為了驗(yàn)證該模型的計(jì)算精度，對(duì)測(cè)試樣本進(jìn)行了模擬計(jì)算，并計(jì)算了模擬結(jié)果與克里金模型預(yù)測(cè)結(jié)果的相對(duì)誤差，其結(jié)果如表3所示。

由表3的計(jì)算結(jié)果可知，通過(guò)對(duì)聚能射流速度模擬結(jié)果與克里金模型預(yù)測(cè)結(jié)果進(jìn)行相對(duì)誤差分析，其最大誤差不超過(guò)2%，可知克里金模型預(yù)測(cè)值與模擬計(jì)算結(jié)果吻合較好，表明克里金模型預(yù)測(cè)聚能射流速度時(shí)具有較高的準(zhǔn)確度。在實(shí)際計(jì)算中，只需要11組樣本的數(shù)值計(jì)算結(jié)果，并利用克里金模型，即可獲得測(cè)試樣本的射流速度，即通過(guò)11次數(shù)值模擬獲得了27組樣本的射流速度結(jié)果，其計(jì)算效率大幅提高。需要說(shuō)明的是，射流長(zhǎng)度等可以用單個(gè)參數(shù)表征的射流性能也可通過(guò)本研究的方法類似計(jì)算，這里不再贅述。

3 結(jié) 論

(1)所建立的自適應(yīng)克里金模型可用于快速預(yù)測(cè)聚能裝藥射流速度，與有限元模擬結(jié)果相比，兩者誤差小于2%，表明克里金模型對(duì)射流速度計(jì)算具有較高的準(zhǔn)確度。

(2)克里金模型能夠用少量的計(jì)算成本，獲得聚能裝藥射流速度，使得計(jì)算效率大幅提升，該方法提供了一條有效的技術(shù)途徑，可用于指導(dǎo)聚能裝藥研究和理論設(shè)計(jì)。