侵徹條件下巖石類材料RHT模型參數敏感性分析

聶錚玥， 彭 永， 陳 榮， 蔡鎮清， 李付剛

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

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

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

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

1 RHT模型簡介及參數分類

1.1 RHT模型介紹

RHT模型由Riedel等[14]提出，可以反映材料的應變率效應、失效面及損傷演化等特點。RHT模型包括狀態方程部分及本構方程部分，如圖1所示。

圖1 RHT模型原理圖

狀態方程采用Herrmann[15]提出的P-α狀態方程，以描述多孔及疏松介質中孔隙從壓碎到壓實的過程。本構模型由3個發展階段組成，隨著應力增加，材料首先經過彈性階段到達彈性屈服面，此后材料發生塑性變形并進入線性強化階段，材料表現出應變硬化特征，直至到達失效面。當等效應力強度超出失效應力強度后，材料開始出現累積損傷量，進入損傷軟化階段，最后到達殘余強度面[14-16]。

根據RHT模型原理，基于參數物理意義，可對34個參數進行初步歸類，如表1所示。

表1 RHT模型參數分類

1.2 參數分類及范圍選取

根據Rankine-Hugoniot方程及Mie-Grüneisen狀態方程可以求得部分參數間關系為

(1)

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

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

(2)

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

表2 主定參數、從屬參數、定值參數分類

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

表3 主定參數參考取值范圍

2 敏感性分析數值實驗設計

2.1 參數矩陣設計

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

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

Smin,1=M+1

(3)

(4)

式中：Smin為最小采樣點數；M為試驗設計的主定參數個數。

因此，考慮到最小采樣點數要求及計算成本，基于上節中的主定參數范圍，利用Opt LHD算法設置了400個采樣點，獲得了圖2所示的參數設計矩陣。

圖2 主定參數設計矩陣

2.2 數值實驗設計

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

圖3 模型尺寸及網格設計(mm)

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

表4 彈體Johnson-cook模型參數[25]

3 RHT模型參數敏感性分析方法

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

3.1 交互效應、主效應分析原理

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

(5)

式中：xixj為參數交互項；εi為誤差項。

對式(5)微分可得

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

(6)

當式(6)中交互項為零或值很小時，可得到各參數的主效應。線性項的主效應為

Mxi=cidxi

(7)

二階項的主效應為

(8)

當交互項不可忽略時，則參數間的交互效應如式(9)所示

Mxixj=cmd(xixj)

(9)

3.2 參數貢獻率計算原理

為了計算各個參數的貢獻率，需建立樣本點和評價指標之間的近似模型。首先通過式(10)將各主定參數值歸一化至[-1,+1]，以消除21個主定參數自身數值大小的影響。

(10)

主定參數的貢獻率采用響應面模型方法來計算，模型可分為多元一次回歸模型、多元二次回歸模型等，表達式分別為

(11)

(12)

(13)

3.3 評價指標選取

如圖4所示，不同參數條件下，彈體侵徹深度、彈體磨蝕程度及靶體損傷程度均不同。因此，分別選擇侵徹深度H、彈體質量損失率及靶體損傷度D作為評價指標，進行參數敏感性分析。

圖4 不同工況下的仿真結果

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

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

4 RHT模型參數敏感性結果

數值模擬獲得的有效工況共計355個。剩余45個工況則因為參數取值不合理，導致計算時靶體單元體中出現復雜聲速，或由于靶體強度過低，彈體將貫穿靶體而不采用。提取各工況的侵徹深度、彈體質量損失率及靶體損傷度作為敏感性分析的評價指標進行后續分析。

4.1 參數交互效應分析

圖6 對于不同評價指標的部分參數間的交互效應圖

因此，對于同一評價指標，不同參數之間的交互效應不同。同時，對于不同評價指標，某兩個參數之間的交互效應也可能存在差異。為了進一步描述參數之間的交互效應，取21個主定參數進行總體參數間交互效應分析。

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

4.2 參數主效應分析

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

圖7 侵徹深度主效應圖

圖8 彈體質量損失率主效應圖

圖9 靶體損傷度主效應圖

4.3 近似模型及參數貢獻率分析

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

圖10 各階響應面模型誤差分析圖

表5 各階響應面模型可決系數R2結果

基于二階響應面模型，分別對侵徹深度、彈體質量損失率及靶體損傷度計算各參數項的貢獻率，并提取絕對值較大的前21項，如圖11所示。結果表明，對于侵徹深度及靶體損傷度而言，參數二次項對評價指標的影響普遍較一次項及交互項大。對于彈體質量損失率而言，貢獻率較大的參數項同樣以二次項居多。

為了較全面地衡量各主定參數的貢獻率情況，按下式計算各個主定參數總體貢獻率

(14)

表6 各主定參數總體貢獻率

5 結 論

本文針對彈體侵徹半無限厚巖石靶體問題，利用優化拉丁超立方設計算法設計了參數試驗矩陣，并進行數值模擬。根據355個有效結果，以侵徹深度、彈體質量損失率及靶體損傷度為評價指標，對21個主定參數進行了參數敏感性分析，并得到了以下結論：

(1) 主定參數間存在一定交互效應，且不同參數之間的交互效應不同。同時，對于不同評價指標，參數間的交互效應亦存在差異。取21個主定參數對侵徹深度、彈體質量損失率及靶體損傷度進行總體參數組合間交互效應分析，其結果均為交互效應不顯著。

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

(3) 對于不同評價指標，各參數的主效應不一致。對于侵徹深度及彈體質量損失率而言，殘余應力強度指數的主效應最顯著。此外，當參數水平變化時，各主定參數間的主效應相對顯著情況將發生變化。