混凝土侵徹數值模擬的影響因素＊

林華令，丁育青，湯文輝

（1.第二炮兵裝備研究院第二研究所，北京100085；2.國防科學技術大學理學院工程物理研究所，湖南 長沙410073）

近年來，有限元數值計算方法在混凝土侵徹問題的分析中得到了廣泛應用，而計算結果的準確性與混凝土材料動態本構模型及參數、網格單元尺寸等因素直接相關。當使用拉格朗日網格進行數值模擬時，由于在侵徹過程中靶板網格會產生畸變，因而，網格銷蝕判據對計算結果也有顯著影響。為考察上述因素對計算結果的影響規律，可以對混凝土侵徹實驗過程進行數值模擬，并將計算結果與實驗結果進行比較，從而為數值模擬工作提供參考，并有效提高計算的準確度。門建兵等［1］針對網格尺寸對混凝土侵徹過程的影響進行了數值計算，結果表明，侵徹彈丸半徑與靶板網格邊長的比值在6.0左右，計算結果較理想。J.Lepp?en［2］針對影響侵徹深度和靶板前表面損傷面積計算結果的相關參數進行了研究，結果表明，網格尺寸和網格銷蝕判據對侵徹深度計算結果影響明顯，但對靶板前表面損傷面積計算結果影響不大，而混凝土損傷計算結果與其動態本構模型中的拉伸描述直接相關。AUTODYN軟件［3］中，混凝土RHT模型提供2種拉伸失效模式，而選擇不同拉伸失效模式對計算結果的影響尚未見報道。本文中，利用AUTODYN軟件，采用混凝土RHT模型，對混凝土侵徹過程進行數值模擬。針對混凝土拉伸失效模式、混凝土靶板網格劃分以及網格銷蝕判據等因素進行考察，并對計算結果進行分析，以期得到以上探討因素對計算結果的影響規律。

1 數值模擬

依據H.Hansson［4］的實驗結果，對混凝土侵徹過程進行數值模擬。該實驗中侵徹子彈為卵形彈，彈徑比為3.0，直徑為75mm，彈體長為225mm，密度為7 830kg／m3，彈體質量約為6.28kg。彈體材料為4340鋼，體積模量為159GPa，剪切模量為81.8GPa，屈服強度為792MPa，入射速度為485m／s。混凝土靶板為圓柱形，長度為2.0m，直徑為1.6m，邊長為150mm的標準立方塊試樣的抗壓強度為40MPa，單軸抗拉強度為2.64MPa，斷裂能為100J／m2。2次重復實驗侵徹深度分別為655、660mm，混凝土靶板表面成坑直徑約為800mm。

考慮到軸對稱性，利用AUTODYN軟件進行二維數值模擬，對彈、靶模型均采用拉格朗日網格進行描述。為了研究靶板網格尺寸對計算結果的影響，將網格形狀取為正方形，劃分的網格尺寸分別為10.0、8.0、5.0和2.5mm，對應的單元數量分別為200×80、250×100、400×160和800×320，網格尺寸取5.0mm時的彈、靶局部網格模型如圖1所示。

圖1 彈靶網格劃分模型Fig.1 Numerical meshes for projectile and target

計算過程中，對子彈采用AUTODYN軟件材料庫中的STEEL 4340模型［3］；對混凝土靶板采用RHT本構模型［5］，該模型綜合考慮了混凝土失效面的壓力相關性、壓縮損傷軟化、應變率效應等特點，同時引入了偏應力張量第三不變量對失效面形狀的影響，并考慮了拉、壓應變率效應的差異。強度描述方面，RHT模型中引入了3個失效面，即最大失效面、彈性屈服失效面和殘余失效面，分別描述混凝土的失效強度、初始屈服強度及殘余強度的變化。RHT模型中采用考慮多孔度的p－α狀態方程［6］。對于混凝土材料，AUTODYN軟件提供了2種拉伸失效模式：靜水壓拉伸失效和主應力拉伸失效。其中，靜水壓拉伸失效判據是當靜水壓達到預設值后判定失效，主應力拉伸失效判據是當某一方向主應力的最大拉伸應力達到預設值后判定失效。由于RHT模型在描述混凝土拉伸失效時存在不足［7］，本文中采用2種拉伸失效模式分別進行計算，進而對2種模式進行比對。計算過程中，混凝土靶板的網格會產生畸變，因此需要采用網格銷蝕判據。AUTODYN軟件提供了多種侵蝕／銷蝕模型，其中，瞬時幾何應變在侵徹計算中較常用，銷蝕判據與混凝土靶板的網格大小直接相關。其次，由于RHT模型關于混凝土強度的描述考慮了應變率，而應變率的計算與網格大小也存在關聯，因此需要對網格劃分以及銷蝕應變對于混凝土侵徹數值計算結果的影響進行考察。

計算過程中，不考慮彈體和混凝土靶板之間的摩擦，考慮靶板混凝土材料的裂紋軟化。根據實驗靶板的混凝土材料參數對部分RHT模型參數進行調整：剪切模量G＝14.33GPa，單軸抗壓強度fc＝33.8MPa，單軸拉／壓強度比ft／fc＝0.078，由文獻［8］計算得出；失效面參數A＝1.929，失效面參數N＝0.764，拉壓子午比Q0＝0.69，脆性韌性轉變參數BQ＝0.004 8，取自文獻［9］；拉伸應變率指數δ＝0.025，殘余強度面參數B＝1.5，殘余強度面參數M＝0.7，取自文獻［10］。其余參數采用AUTODYN軟件材料庫中CONC－35MPa中的默認參數。拉伸失效模式分別選用靜水壓拉伸失效和主應力拉伸失效，主應力拉伸閾值取5MPa［11］，銷蝕應變分別取1.0、1.5、2.0、2.5和3.0。

2 計算結果及分析

2.1 靶板網格劃分對計算結果的影響

為了單獨考察靶板網格劃分對計算結果的影響，采用固定模型參數（銷蝕應變取1.5）、僅將網格尺寸作為變量的方法，針對不同網格劃分的靶板模型分別進行計算，計算結果如圖2所示。

圖2 網格尺寸對侵徹深度的影響Fig.2 Effects of meshing on depth of penetration

采用靜水壓拉伸失效模式時，網格尺寸l對侵徹深度s計算結果的影響如圖2（a）所示，可以看出：侵徹深度隨著網格尺寸增大而減小，單元網格尺寸取10.0mm時，侵徹深度計算結果最小；網格尺寸分別取8.0、5.0和2.5mm時，侵徹深度計算結果差別不大，與實驗結果較接近，偏小約7%；網格尺寸取2.5mm時，由于計算步長減小，計算耗時明顯增加，而增加的計算時間對計算精度沒有明顯改進。采用主應力拉伸失效模式時，網格尺寸對侵徹深度計算結果的影響如圖2（b）所示，可以看出：侵深計算結果與網格尺寸沒有明顯的線性關系；單元網格尺寸取2.5mm時，侵徹深度計算結果與實驗結果差別最大；單元網格尺寸取5.0mm時，計算結果與實驗結果較接近，偏小約2%。

圖3 網格尺寸對靶板損傷區域的影響Fig.3 Effects of meshing on crater size

網格劃分對靶板表面損傷計算結果的影響如圖3所示，圖中箭頭指向位置分別對應損傷半徑與侵徹深度，將靶板前表面損傷區域邊界定為表面成坑半徑的計算結果［10］。由于單元網格尺寸取8.0mm 與取5.0mm 時的計算結果較接近，圖中未列出網格尺寸取8.0mm時的計算結果。從圖3可以看出：采用靜水壓拉伸失效模式，網格尺寸分別取2.5和5.0mm時，損傷計算結果差異不大，且與實驗結果接近，而網格尺寸取10.0mm時，損傷計算結果明顯變小；采用主應力拉伸失效模式時，網格尺寸對損傷計算結果影響不明顯，3種網格的計算結果比實驗結果偏小約25%；對于2種拉伸失效模式，靶板網格尺寸取5.0mm時，計算結果均與實驗結果最接近，且效率較高。

2.2 銷蝕應變對計算結果的影響

為了單獨考察銷蝕應變對計算結果的影響，固定靶板單元網格尺寸為5.0mm，選用瞬時銷蝕失效模型，僅將銷蝕應變作為唯一的變量，計算結果如圖4所示。

圖4 銷蝕應變對侵徹深度的影響Fig.4 Effects of erosion strain on depth of penetration

采用靜水壓拉伸失效模式時，銷蝕應變對侵徹深度計算結果的影響如圖4（a）所示，可以看出，銷蝕應變對計算結果影響明顯：銷蝕應變越大，侵徹深度越小；銷蝕應變取2.0時，與實驗結果吻合較好；銷蝕應變取2.5和3.0時，計算結果與實驗結果差異較大。采用主應力拉伸失效模式時，銷蝕應變對侵徹深度計算結果的影響如圖4（b）所示，可以看出：侵徹深度與銷蝕應變成非線性變化關系，銷蝕應變取1.5時，計算結果與實驗結果較接近；銷蝕應變取其他值時的計算結果均與實驗結果差異很大，誤差大于30%。

銷蝕應變對靶板表面損傷計算結果的影響如圖5所示，其中，銷蝕應變取1.5時的損傷圖與圖3中單元網格尺寸取5.0mm時的損傷一致。可以看出，損傷面積隨著銷蝕應變的增大而增大。采用靜水壓拉伸失效模式時，銷蝕應變取2.0的損傷計算結果與實驗結果較接近。采用主應力拉伸失效模式時，銷蝕應變取1.0、1.5和2.0時的損傷計算結果之間差異不大，比實驗結果均偏小約20%，而當銷蝕應變取2.5和3.0時，損傷計算結果增大較明顯，但比實驗結果偏大。

圖5 銷蝕應變對損傷區域的影響Fig.5 Effects of erosion strain on crater size

3 結 論

利用AUTODYN軟件，采用混凝土RHT模型，對卵形彈垂直侵徹混凝土靶的過程進行數值模擬。通過設置不同的拉伸失效模式、網格劃分和銷蝕應變，進而考察其對混凝土侵徹數值計算結果的影響。計算結果表明，以上3種因素對計算結果的影響均較顯著。選擇靜水壓失效模式能較好地對混凝土侵徹過程進行數值模擬，但需結合混凝土靶板的實際情況，首先確定較準確的模型參數。選擇主應力拉伸失效模式，侵徹深度和損傷計算結果不能同時較好地模擬實驗結果，這與AUTODYN軟件中裂紋軟化模型默認設置為線性關系有關，而實驗結果表明，用雙線性關系描述混凝土的裂紋軟化更合適［10－12］。就本算例來說，分別選用2種不同的拉伸失效模式，靶板單元網格尺寸取5.0mm、侵蝕應變取1.5時，計算結果與實驗結果均較接近；而當選用靜水壓拉伸失效模式，靶板單元網格尺寸取5.0mm、侵蝕應變取2.0時，計算結果與實驗結果較接近。

［1］門建兵，隋樹元，蔣建偉，等.網格對混凝土侵徹數值模擬的影響［J］.北京理工大學學報，2005，25（8）：659－662.Men Jian－bing，Sui Shu－yuan，Jiang Jian－wei，et al.Mesh dependency for numerical simulation of concrete penetration［J］.Transactions of Beijing Institute of Technology，2005，25（8）：659－662.

［2］Lepp?en J.Dynamic behaviour of concrete structures subjected to blast and fragment impacts［D］.G?teborg，Swe－den：Chalmers University of Technology，2002：42－53.

［3］Century Dynamics Inc.AUTODYN manuals：Version 5［M］.Concord，CA，USA：Century Dynamics Inc，2004.

［4］Hansson H.Numerical simulation of concrete penetration［R］.FOA Report 98－0081－311－SE，1998.

［5］Riedel W，Thoma K，Hiermaier S.Penetration of reinforced concrete by BETA－B－500－numerical analysis using a new macroscopic concrete model for hydrocodes［C］∥Proceedings of 9th International Symposium on Interaction of the Effect of Munitions with Structures.Germany：Berlin－Strausberg，1999：315－322.

［6］Herrmann W.Constitutive equation for the dynamic compaction of ductile porous materials［J］.Journal of Applied Physics，1969，40（6）：2490－2499.

［7］Tu Zheng－guo，Lu Yong.Evaluation of typical concrete material models used in hydrocodes for high dynamic response simulations［J］.International Journal of Impact Engineering，2009，36（1）：132－146.

［8］Comite Euro－International de Beton.CEB－FIP model code 1990［M］.Trowbridge，UK：Redwood Books，1993.

［9］張若棋，丁育青，湯文輝，等.混凝土 HJC、RHT本構模型的失效強度參數［J］.高壓物理學報，2011，25（1）：15－22.Zhang Ruo－qi，Ding Yu－qing，Tang Wen－hui，et al.The failure strength parameters of concrete HJC and RHT constitutive model［J］.Chinese Journal of High Pressure Physics，2011，25（1）：15－22.

［10］Nystr?m U，Gylltoft K.Comparative numerical studies of projectile impacts on plain and steel－fibre reinforced concrete［J］.International Journal of Impact Engineering，2011，38（2／3）：95－105.

［11］Tu Zheng－guo，Lu Yong.Modifications of RHT material model for improved numerical simulation of dynamic response of concrete［J］.International Journal of Impact Engineering，2010，37（10）：1072－1082.

［12］Lepp?en J.Concrete subjected to projectile and fragment impacts：Modelling of crack softening and strain rate dependency in tension［J］.International Journal of Impact Engineering，2006，32（11）：1828－1841.