混凝土侵彻数值模拟的影响因素＊

林华令，丁育青，汤文辉

（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.