聶錚玥， 彭 永， 陳 榮， 蔡鎮(zhèn)清， 李付剛

(國防科技大學(xué) 文理學(xué)院，長沙 410073)

RHT(Riedel-Hiermaier-Thoma)模型廣泛應(yīng)用于巖石、混凝土類材料在高應(yīng)變率、大變形問題中的動(dòng)態(tài)力學(xué)響應(yīng)模擬[1]，國內(nèi)外學(xué)者基于RHT模型開展的爆炸沖擊、彈體侵徹等研究都取得了比較理想的結(jié)果[2-4]。目前多數(shù)研究的重點(diǎn)在于，針對(duì)混凝土材料，基于實(shí)驗(yàn)數(shù)據(jù)對(duì)部分參數(shù)或模型進(jìn)行修正與改進(jìn)[5-8]。然而，當(dāng)RHT模型應(yīng)用于巖石材料時(shí)，由于不同種類巖石間的力學(xué)性能差異明顯，一般需通過理論分析與力學(xué)實(shí)驗(yàn)相結(jié)合的方法來確定模型參數(shù)，不宜直接采用某種特定巖石或混凝土材料的參考值[9]。此外，RHT模型參數(shù)多達(dá)34個(gè)，獲得參數(shù)所需力學(xué)實(shí)驗(yàn)包括物理性質(zhì)測(cè)定、靜態(tài)力學(xué)及動(dòng)態(tài)力學(xué)等多種類型，使得參數(shù)確定比較復(fù)雜。因此，對(duì)模型中參數(shù)進(jìn)行敏感性分析，甄別敏感性較強(qiáng)的參數(shù)，可以優(yōu)化模型參數(shù)確定方案，提高數(shù)值模擬的效率。

對(duì)于RHT模型參數(shù)敏感性的研究相對(duì)較少，李洪超等[10-11]采用LS-DYNA數(shù)值軟件并通過正交試驗(yàn)的方法，模擬分析了RHT模型中的部分復(fù)雜參數(shù)對(duì)大理巖及紅砂巖霍普金森桿(split Hopkinson pressure bar, SHPB)壓縮實(shí)驗(yàn)應(yīng)力-應(yīng)變曲線結(jié)果的敏感性，通過極差分析，給出了影響明顯的材料模型參數(shù)排序。該方法實(shí)驗(yàn)設(shè)計(jì)合理，且能夠在全局上獲得擬合SHPB沖擊曲線的最佳解。但由于SHPB沖擊曲線結(jié)果將受材料物理性質(zhì)影響而具有一定離散性，且不同應(yīng)變率條件下獲得的SHPB沖擊曲線也將不同，因此該方法有一定局限性。劉殿柱等[12]同樣利用正交試驗(yàn)法，以侵徹深度為評(píng)價(jià)指標(biāo)，利用Autodyn軟件對(duì)彈丸侵徹混凝土問題進(jìn)行了部分參數(shù)的敏感性分析。該方法簡(jiǎn)明直觀，但根據(jù)侵徹深度變化趨勢(shì)圖對(duì)敏感性進(jìn)行排序時(shí)，由于部分參數(shù)影響相對(duì)不顯著，因此結(jié)果具有一定主觀性。辛健[13]通過模擬TNT在土中的爆炸過程，采用單因素變量法，得到了以損傷面積為評(píng)價(jià)指標(biāo)的參數(shù)敏感度結(jié)果。該方法簡(jiǎn)潔直觀，但參數(shù)水平設(shè)置情況通常有限，且RHT模型中參數(shù)并非完全獨(dú)立，因此分析結(jié)果可能存在一定誤差。為客觀評(píng)價(jià)模型參數(shù)的敏感性差異，需要對(duì)計(jì)算參數(shù)矩陣進(jìn)行合理設(shè)計(jì)與分析。同時(shí)，上述研究結(jié)果表明，針對(duì)TNT爆炸、SHPB試驗(yàn)與侵徹問題，RHT模型材料參數(shù)的敏感性不同。因此，針對(duì)特定問題，需要進(jìn)行相應(yīng)的參數(shù)敏感性分析，不能直接采用其他物理問題的敏感性分析結(jié)果。

針對(duì)實(shí)心彈體中低速正侵徹單層半無限厚巖石靶體問題，本研究基于ISIGHT軟件中的優(yōu)化拉丁超立方設(shè)計(jì) (optimal Latin hypercube design，Opt LHD)算法設(shè)計(jì)了參數(shù)試驗(yàn)矩陣，保證了設(shè)計(jì)矩陣具有較優(yōu)的空間填充性。采用LS-DYNA軟件對(duì)侵徹過程進(jìn)行了大量數(shù)值模擬，獲得了不同參數(shù)組合條件下的仿真結(jié)果。然后以彈體侵徹深度、彈體質(zhì)量損失率及靶體損傷度為評(píng)價(jià)指標(biāo)，分析了模型參數(shù)的交互效應(yīng)及主效應(yīng)。最后通過建立侵徹深度、彈體質(zhì)量損失率及靶體損傷度與參數(shù)間的響應(yīng)面模型，定量計(jì)算了RHT模型中參數(shù)的貢獻(xiàn)率。

1 RHT模型簡(jiǎn)介及參數(shù)分類

1.1 RHT模型介紹

RHT模型由Riedel等[14]提出，可以反映材料的應(yīng)變率效應(yīng)、失效面及損傷演化等特點(diǎn)。RHT模型包括狀態(tài)方程部分及本構(gòu)方程部分，如圖1所示。

圖1 RHT模型原理圖

狀態(tài)方程采用Herrmann[15]提出的P-α狀態(tài)方程，以描述多孔及疏松介質(zhì)中孔隙從壓碎到壓實(shí)的過程。本構(gòu)模型由3個(gè)發(fā)展階段組成，隨著應(yīng)力增加，材料首先經(jīng)過彈性階段到達(dá)彈性屈服面，此后材料發(fā)生塑性變形并進(jìn)入線性強(qiáng)化階段，材料表現(xiàn)出應(yīng)變硬化特征，直至到達(dá)失效面。當(dāng)?shù)刃?yīng)力強(qiáng)度超出失效應(yīng)力強(qiáng)度后，材料開始出現(xiàn)累積損傷量，進(jìn)入損傷軟化階段，最后到達(dá)殘余強(qiáng)度面[14-16]。

根據(jù)RHT模型原理，基于參數(shù)物理意義，可對(duì)34個(gè)參數(shù)進(jìn)行初步歸類，如表1所示。

表1 RHT模型參數(shù)分類

1.2 參數(shù)分類及范圍選取

根據(jù)Rankine-Hugoniot方程及Mie-Grüneisen狀態(tài)方程可以求得部分參數(shù)間關(guān)系為

(1)

式中：c0為材料聲速；s為材料參數(shù)，其值可以通過平板撞擊試驗(yàn)中測(cè)定的材料沖擊波速度和波后粒子速度結(jié)果線性擬合確定。其余參數(shù)間亦存在下列關(guān)系[17-18]

pel=fc/3
βc=4/(20+3fc)
βt=2/(20+fc)

(2)

因此，可將c0，s，fc作為主定參數(shù)，而A1，A2，A3，B0，B1，T1，pel，βc，βt作為其從屬參數(shù)處理。綜合各參數(shù)定義及式(1)、式(2)，主定參數(shù)、從屬參數(shù)及定值參數(shù)的具體分類如表2所示，其中主定參數(shù)共計(jì)21個(gè)。

表2 主定參數(shù)、從屬參數(shù)、定值參數(shù)分類

本研究針對(duì)巖石類材料，參考花崗巖、玄武巖、砂巖、千枚巖、大理巖、片麻巖共六種巖石的相關(guān)指標(biāo)值[19]，確定了部分主定參數(shù)的取值范圍。對(duì)于材料參數(shù)s的取值范圍則根據(jù)巖石材料的平板撞擊試驗(yàn)結(jié)果確定[20-23]。對(duì)于無法直接確定的其余主定參數(shù)，則利用C30/37混凝土RHT模型參數(shù)值作為基準(zhǔn)值，并設(shè)置參數(shù)變化范圍為基準(zhǔn)值的-80%～80%。所有主定參數(shù)的取值范圍如表3所示。

表3 主定參數(shù)參考取值范圍

2 敏感性分析數(shù)值實(shí)驗(yàn)設(shè)計(jì)

2.1 參數(shù)矩陣設(shè)計(jì)

ISIGHT是一款可以靈活設(shè)計(jì)并調(diào)用各項(xiàng)軟件進(jìn)行計(jì)算，進(jìn)而對(duì)結(jié)果進(jìn)行后處理的計(jì)算機(jī)輔助優(yōu)化軟件[24]。本研究利用ISIGHT中的試驗(yàn)設(shè)計(jì)(design of experiments, DOE)功能，采用優(yōu)化拉丁超立方設(shè)計(jì)算法(Opt LHD)進(jìn)行參數(shù)矩陣設(shè)計(jì)，從而保證試驗(yàn)點(diǎn)具有較好的空間填充性及均勻度。

根據(jù)RHT模型中21個(gè)主定參數(shù)設(shè)計(jì)參數(shù)矩陣，為滿足后續(xù)的交互效應(yīng)及主效應(yīng)分析，以及建立響應(yīng)面模型的精度需求，參數(shù)水平個(gè)數(shù)的設(shè)計(jì)需滿足最小采樣點(diǎn)數(shù)。對(duì)于一階、二階響應(yīng)面模型，其最小采樣點(diǎn)數(shù)分別通過式(3)、式(4)計(jì)算

Smin,1=M+1

(3)

(4)

式中：Smin為最小采樣點(diǎn)數(shù)；M為試驗(yàn)設(shè)計(jì)的主定參數(shù)個(gè)數(shù)。

因此，考慮到最小采樣點(diǎn)數(shù)要求及計(jì)算成本，基于上節(jié)中的主定參數(shù)范圍，利用Opt LHD算法設(shè)置了400個(gè)采樣點(diǎn)，獲得了圖2所示的參數(shù)設(shè)計(jì)矩陣。

圖2 主定參數(shù)設(shè)計(jì)矩陣

2.2 數(shù)值實(shí)驗(yàn)設(shè)計(jì)

利用LS-DYNA有限元軟件對(duì)實(shí)心彈體中低速正侵徹半無限厚單層巖石靶體問題進(jìn)行仿真。為節(jié)省計(jì)算時(shí)間，采用1/4對(duì)稱模型進(jìn)行仿真實(shí)驗(yàn)，尺寸及網(wǎng)格設(shè)計(jì)如圖3所示，彈體在靶體中心處以600 m/s的初速度垂直向下侵徹。靶體的上表面設(shè)置為自由邊界，下表面設(shè)置為無反射邊界，以近似模擬半無限厚條件。

圖3 模型尺寸及網(wǎng)格設(shè)計(jì)(mm)

彈體的材料選擇為硬度較大的 30CrMnSiNi2A 鋼，使用Johnson-Cook模型描述其材料力學(xué)性能，模型參數(shù)如表4所示。靶體材料采用RHT模型，其中主定參數(shù)為參數(shù)設(shè)計(jì)矩陣值，其余參數(shù)確定如1.2節(jié)所述。

表4 彈體Johnson-cook模型參數(shù)[25]

3 RHT模型參數(shù)敏感性分析方法

為了定量描述參數(shù)敏感性大小以及參數(shù)變化對(duì)評(píng)價(jià)指標(biāo)(即響應(yīng)量)的影響規(guī)律，需對(duì)參數(shù)的交互效應(yīng)、主效應(yīng)以及參數(shù)對(duì)響應(yīng)量的貢獻(xiàn)率進(jìn)行分析。參數(shù)的交互效應(yīng)可以反映某一參數(shù)變化對(duì)響應(yīng)量的貢獻(xiàn)隨其他參數(shù)的取值不同而發(fā)生變化的情況。參數(shù)交互效應(yīng)顯著則此時(shí)參數(shù)作用非獨(dú)立。當(dāng)交互效應(yīng)不顯著時(shí)，則通過參數(shù)主效應(yīng)分析，可得到參數(shù)對(duì)響應(yīng)量的影響規(guī)律。通過建立參數(shù)與響應(yīng)量之間的近似模型，進(jìn)一步計(jì)算各參數(shù)貢獻(xiàn)率，評(píng)價(jià)參數(shù)的敏感性大小。

3.1 交互效應(yīng)、主效應(yīng)分析原理

參數(shù)交互效應(yīng)及主效應(yīng)的計(jì)算是先通過式(5)得到主定參數(shù)xi與響應(yīng)量y之間多項(xiàng)式(以二次多項(xiàng)式為例)

(5)

式中：xixj為參數(shù)交互項(xiàng)；εi為誤差項(xiàng)。

對(duì)式(5)微分可得

dy=c1dx1+…+2cn+1x1dx1+…+
c2n+1d(x1x2)+…+εi

(6)

當(dāng)式(6)中交互項(xiàng)為零或值很小時(shí)，可得到各參數(shù)的主效應(yīng)。線性項(xiàng)的主效應(yīng)為

Mxi=cidxi

(7)

二階項(xiàng)的主效應(yīng)為

(8)

當(dāng)交互項(xiàng)不可忽略時(shí)，則參數(shù)間的交互效應(yīng)如式(9)所示

Mxixj=cmd(xixj)

(9)

3.2 參數(shù)貢獻(xiàn)率計(jì)算原理

為了計(jì)算各個(gè)參數(shù)的貢獻(xiàn)率，需建立樣本點(diǎn)和評(píng)價(jià)指標(biāo)之間的近似模型。首先通過式(10)將各主定參數(shù)值歸一化至[-1,+1]，以消除21個(gè)主定參數(shù)自身數(shù)值大小的影響。

(10)

主定參數(shù)的貢獻(xiàn)率采用響應(yīng)面模型方法來計(jì)算，模型可分為多元一次回歸模型、多元二次回歸模型等，表達(dá)式分別為

(11)

(12)

(13)

3.3 評(píng)價(jià)指標(biāo)選取

如圖4所示，不同參數(shù)條件下，彈體侵徹深度、彈體磨蝕程度及靶體損傷程度均不同。因此，分別選擇侵徹深度H、彈體質(zhì)量損失率及靶體損傷度D作為評(píng)價(jià)指標(biāo)，進(jìn)行參數(shù)敏感性分析。

圖4 不同工況下的仿真結(jié)果

其中，侵徹深度根據(jù)計(jì)算結(jié)果直接獲得。彈體的質(zhì)量損失率通過提取彈體單元體質(zhì)量變化曲線計(jì)算得到。而對(duì)于靶體損傷度的表征，則是通過將侵徹過程中已經(jīng)刪除的單元損傷度定義為1，并提取剩余所有單元的損傷度進(jìn)行累加獲得。典型靶體的全模型損傷分布如圖5所示。

圖5 工況41全模型的損傷情況

4 RHT模型參數(shù)敏感性結(jié)果

數(shù)值模擬獲得的有效工況共計(jì)355個(gè)。剩余45個(gè)工況則因?yàn)閰?shù)取值不合理，導(dǎo)致計(jì)算時(shí)靶體單元體中出現(xiàn)復(fù)雜聲速，或由于靶體強(qiáng)度過低，彈體將貫穿靶體而不采用。提取各工況的侵徹深度、彈體質(zhì)量損失率及靶體損傷度作為敏感性分析的評(píng)價(jià)指標(biāo)進(jìn)行后續(xù)分析。

4.1 參數(shù)交互效應(yīng)分析

圖6 對(duì)于不同評(píng)價(jià)指標(biāo)的部分參數(shù)間的交互效應(yīng)圖

因此，對(duì)于同一評(píng)價(jià)指標(biāo)，不同參數(shù)之間的交互效應(yīng)不同。同時(shí)，對(duì)于不同評(píng)價(jià)指標(biāo)，某兩個(gè)參數(shù)之間的交互效應(yīng)也可能存在差異。為了進(jìn)一步描述參數(shù)之間的交互效應(yīng)，取21個(gè)主定參數(shù)進(jìn)行總體參數(shù)間交互效應(yīng)分析。

SPSS(statistical product and service solutions)軟件是一個(gè)功能強(qiáng)大的統(tǒng)計(jì)分析軟件，通過SPSS軟件可以對(duì)參數(shù)交互效應(yīng)進(jìn)行計(jì)算，并獲得相應(yīng)的顯著性結(jié)果(significance，Sig)，當(dāng)Sig值小于0.05時(shí)，說明此時(shí)參數(shù)交互效應(yīng)顯著，參數(shù)間將互相影響，反之則不顯著。結(jié)果表明：當(dāng)評(píng)價(jià)指標(biāo)為侵徹深度時(shí)，21個(gè)主定參數(shù)組合的交叉效應(yīng)Sig值為0.296>0.05；當(dāng)評(píng)價(jià)指標(biāo)為彈體質(zhì)量損失率時(shí)，交叉效應(yīng)Sig值為0.397>0.05；當(dāng)評(píng)價(jià)指標(biāo)為靶體損傷度時(shí)，交叉效應(yīng)Sig值為0.995>0.05。因此，對(duì)于侵徹深度、彈體質(zhì)量損失率與靶體損傷度而言，總體主定參數(shù)組合間的交互效應(yīng)均不顯著。

4.2 參數(shù)主效應(yīng)分析

為了研究參數(shù)變化對(duì)評(píng)價(jià)指標(biāo)的影響規(guī)律，針對(duì)各評(píng)價(jià)指標(biāo)進(jìn)行參數(shù)主效應(yīng)分析，結(jié)果如圖7～圖9所示。其中，橫坐標(biāo)表示參數(shù)取值從低水平(對(duì)應(yīng)1)依次變化至高水平(對(duì)應(yīng)2)，縱坐標(biāo)表示計(jì)算得到的評(píng)價(jià)指標(biāo)的統(tǒng)計(jì)值。對(duì)于侵徹深度與靶體損傷度，各主定參數(shù)的主效應(yīng)均呈非線性及非單調(diào)性變化。當(dāng)各主定參數(shù)從低水平值變化至高水平值時(shí)，侵徹深度及靶體損傷度均將先增大后減小。對(duì)于彈體質(zhì)量損失率，除了材料參數(shù)s外，其余參數(shù)對(duì)指標(biāo)的主效應(yīng)亦呈非線性及非單調(diào)性變化。而當(dāng)s值增加時(shí)，彈體的質(zhì)量統(tǒng)計(jì)損失率將降低。

圖7 侵徹深度主效應(yīng)圖

圖8 彈體質(zhì)量損失率主效應(yīng)圖

圖9 靶體損傷度主效應(yīng)圖

4.3 近似模型及參數(shù)貢獻(xiàn)率分析

分別建立一階、二階、三階響應(yīng)面模型，并利用355個(gè)侵徹深度、彈體質(zhì)量損失率及靶體損傷度結(jié)果進(jìn)行誤差分析，結(jié)果如圖10所示。其中，平行于橫坐標(biāo)的水平線表示用于計(jì)算的所有樣本點(diǎn)的評(píng)價(jià)指標(biāo)的實(shí)際平均值，對(duì)角線表示實(shí)際值與模型預(yù)測(cè)值相等的情況，點(diǎn)越靠近對(duì)角線則表明模型預(yù)測(cè)越精確。同時(shí)，為了定量評(píng)價(jià)各響應(yīng)面模型對(duì)所有實(shí)際評(píng)價(jià)指標(biāo)值的回歸擬合效果，計(jì)算各模型的可決系數(shù)R2(coefficient of determination)如表5所示，當(dāng)R2越接近1，則回歸擬合效果越好。由圖10可知，對(duì)于侵徹深度、彈體質(zhì)量損失率及靶體損傷度，均為一階響應(yīng)面模型誤差顯著大于二階與三階響應(yīng)面模型，即主定參數(shù)與評(píng)價(jià)指標(biāo)間的函數(shù)關(guān)系具有較強(qiáng)的非線性。同時(shí)，由于二階模型與三階模型的可決系數(shù)差異較小，考慮到三階響應(yīng)面模型較二階模型更為復(fù)雜，故選擇二階響應(yīng)面模型進(jìn)行后續(xù)參數(shù)貢獻(xiàn)率計(jì)算。

圖10 各階響應(yīng)面模型誤差分析圖

表5 各階響應(yīng)面模型可決系數(shù)R2結(jié)果

基于二階響應(yīng)面模型，分別對(duì)侵徹深度、彈體質(zhì)量損失率及靶體損傷度計(jì)算各參數(shù)項(xiàng)的貢獻(xiàn)率，并提取絕對(duì)值較大的前21項(xiàng)，如圖11所示。結(jié)果表明，對(duì)于侵徹深度及靶體損傷度而言，參數(shù)二次項(xiàng)對(duì)評(píng)價(jià)指標(biāo)的影響普遍較一次項(xiàng)及交互項(xiàng)大。對(duì)于彈體質(zhì)量損失率而言，貢獻(xiàn)率較大的參數(shù)項(xiàng)同樣以二次項(xiàng)居多。

為了較全面地衡量各主定參數(shù)的貢獻(xiàn)率情況，按下式計(jì)算各個(gè)主定參數(shù)總體貢獻(xiàn)率

(14)

表6 各主定參數(shù)總體貢獻(xiàn)率

5 結(jié) 論

本文針對(duì)彈體侵徹半無限厚巖石靶體問題，利用優(yōu)化拉丁超立方設(shè)計(jì)算法設(shè)計(jì)了參數(shù)試驗(yàn)矩陣，并進(jìn)行數(shù)值模擬。根據(jù)355個(gè)有效結(jié)果，以侵徹深度、彈體質(zhì)量損失率及靶體損傷度為評(píng)價(jià)指標(biāo)，對(duì)21個(gè)主定參數(shù)進(jìn)行了參數(shù)敏感性分析，并得到了以下結(jié)論：

(1) 主定參數(shù)間存在一定交互效應(yīng)，且不同參數(shù)之間的交互效應(yīng)不同。同時(shí)，對(duì)于不同評(píng)價(jià)指標(biāo)，參數(shù)間的交互效應(yīng)亦存在差異。取21個(gè)主定參數(shù)對(duì)侵徹深度、彈體質(zhì)量損失率及靶體損傷度進(jìn)行總體參數(shù)組合間交互效應(yīng)分析，其結(jié)果均為交互效應(yīng)不顯著。

(2) 對(duì)于侵徹深度與靶體損傷度，各主定參數(shù)的主效應(yīng)曲線均呈非線性及非單調(diào)性變化。隨著參數(shù)水平值增加，侵徹深度及靶體損傷度的統(tǒng)計(jì)值均將先增大后減小。而對(duì)于彈體質(zhì)量損失率，材料參數(shù)s的值增加時(shí)，彈體質(zhì)量統(tǒng)計(jì)損失率將降低。其余參數(shù)的主效應(yīng)則與侵徹深度和靶體損傷度相似。

(3) 對(duì)于不同評(píng)價(jià)指標(biāo)，各參數(shù)的主效應(yīng)不一致。對(duì)于侵徹深度及彈體質(zhì)量損失率而言，殘余應(yīng)力強(qiáng)度指數(shù)的主效應(yīng)最顯著。此外，當(dāng)參數(shù)水平變化時(shí)，各主定參數(shù)間的主效應(yīng)相對(duì)顯著情況將發(fā)生變化。