基于SPH-FEM方法的地下结构侵彻爆炸数值模拟研究

2023-04-20张文堂孙惠香袁英杰孙慧颖

弹箭与制导学报 2023年1期

张文堂，孙惠香，袁英杰，孙慧颖，康婷

(1 空军工程大学航空工程学院，陕西西安 710038；2 95445部队，云南大理 672100；3 中天西北建设集团有限公司，陕西西安 710077)

0 引言

钻地武器打击地下掩体一般分为侵彻和爆炸两个阶段,是高烈度、大规模、高时效的复杂毁伤过程。侵爆试验受成本高、难度大、危险性大和测量困难限制,难以实时开展并满足研究需求。数值模拟在侵彻、爆炸领域得到越来越广泛的应用。地下结构在单一侵彻或爆炸作用下的响应分析已有大量研究[1-8]。随着Swegle等[9]将SPH方法应用到水下爆炸;Libersky等[10]将SPH方法应用到材料冲击试验模拟开始,SPH方法解决了许多固体侵彻、爆炸问题。强洪夫等[11]将SPH方法应用到低速碰撞花岗岩领域并进行了仿真分析,通过对比发现SPH方法具有良好的适应性和准确性;王海兵等[12]在试验基础上采用SPH算法较好的模拟了不同弹速撞击花岗岩的试验结果;杨广栋等[13]采用SPH算法先模拟弹体侵彻,再采用Lagrange算法模拟装药爆炸作用下混凝土靶的毁伤效应,提出初始侵彻毁伤对混凝土内爆的重要影响,但侵彻成坑后填充炸药的方法与实际侵爆连续过程仍有差距;胡英国等[14]应用SPH-FEM方法解决了岩体爆破近区精细化数值模拟的问题,与试验监测结果拟合度较好。现有研究只针对单一侵彻或爆炸作用下结构动力响应,或局限于岩石、混凝土靶单一材料的毁伤分析,对大型地下洞室多材料结构在遭受精确制导武器侵彻爆炸连续作用下结构动力响应研究较少。

在SPH算法中,由于质点之间不存在网格关系,因此可避免大变形时网格扭曲造成的精度破坏问题,为解决FEM算法侵爆近区高烈度、大变形以及连续性模拟难题,文中基于SPH-FEM方法建立了钻地弹侵爆地下直墙拱结构的模拟模型,模拟了钻地弹打击地下结构的侵彻爆炸全过程,通过爆腔应力对比了侵彻爆炸与卸荷爆破的差异,分析了直墙拱结构内部应力情况和钻地武器侵彻爆炸联合作用下毁伤模式。

1 SPH-FEM耦合方法

1.1 SPH算法与FEM算法

SPH 算法的基本原理是将连续的流体或固体用相互作用的质点组来描述,各个质点上承载各种物理量,通过求解质点组的动力学方程和跟踪每个质点的运动轨迹,求得整个系统的力学行为。

SPH算法是模拟流体流动的一种拉格朗日型粒子方法,通过任意分布的粒子来求解具有各种边界条件的积分方程或偏微分方程。每个粒子代表一个具有独自物理特性的插值点,用规则的内插函数计算出其力学特性。

(1)

(2)

式中:<>为近似符号;i,j为粒子编号;m,ρ分别为粒子的质量、密度;N为计算区域内粒子总数;Wij=W(|xi-xj|,h)为核函数,h为光滑长度,通常选用三次样条函数。

SPH算法的无网格自适应特性避免了网格畸变计算终止问题,但由于算法耗时巨大,难以应用于大规模工况的模拟。

FEM算法即传统有限元方法,在处理小变形问题时具有较高精度和计算速度。但侵彻爆炸近区具有大变形和非连续特征,采用FEM方法,容易发生网格畸变,导致结果失真或计算终止等问题。

两种算法的优缺点,运用SPH-FEM耦合算法,在侵爆近区使用SPH粒子模拟钻地弹和岩体的变形与破碎,解决网格畸变问题,在侵爆远区采用常规FEM方法保证计算精度、缩短计算时间。

1.2 SPH-FEM耦合方式

钻地弹打击过程分为侵彻和爆炸两个阶段,为实现侵彻和爆炸冲击波产生剥落物飞溅及抛射过程,将弹体和侵彻爆炸近区岩体设置成SPH粒子,以解决FEM因变形过大导致网格畸变运算终止的问题。钻地弹打击地下工程模拟工况体量庞大,考虑计算时间成本和计算精度要求,在地下工程除侵彻爆炸近区外采用有限元网格,提高计算效率。

运用ANSYS/LS-DYNA求解器实现无网格Lagrange求解器和Lagrange处理器的耦合。首先运用LS-Prepost处理器将已经建立好的数值模型中弹体和侵彻爆炸近区岩体的网格单元转换为SPH粒子,经过必要的调试和运行[15],SPH粒子离散间距设置为网格单元尺寸的一半。然后通过设置关键字*CONTACT_TIED_NODES_TO_SURFACE_CONSTRAINED_OFFSET将岩体损伤近区与远区耦合,实现应力波从近区向远区的有效传递,通过比对相邻粒子与单元的应力时程获取最佳参数,实现SPH-FEM耦合计算。

2 计算工况及本构模型

2.1 计算工况

某地下洞库为直墙拱结构,埋深10 m,周围为三级花岗岩,纵深10 m,等间距布置有2 m长加固锚杆,直径2.2 cm;被覆拱为钢筋混凝土结构,净跨为14 m,内拱矢高4.5 m,内墙高2 m,钢筋直径为1.8 cm。

2.2 模型建立

运用传统有限元方法完成建模,模型埋深10 m,纵深10 m,在侵彻爆炸近区4 m×4 m×7 m立方体区域将FEM网格置换为SPH粒子,为避免影响计算精度,对SPH粒子附近及拱结构网格适当加密;为实现地下半无限域模拟,模型四周设置无反射边界条件。具体布置如图1所示。

图1 地下洞库口部截面图(单位:cm)Fig.1 Section of underground cavern mouth(unit:cm)

选择GBU-28钻地弹作为侵彻爆炸弹[16],弹丸总长5.84 m,圆柱部直径37 cm,弹体壁厚12.77 mm,导弹总重2 300 kg,其中装药406 kg,CRH=R/D=3.0。建立1∶1数值模型如图2所示,将模型置换为SPH粒子。

图2 钻地弹模型Fig.2 Model of ground penetrating projectile

2.3 本构模型及参数

钻地弹弹体材料为93钨合金弹,选择Johnson-Cook本构方程描述[17],参数选择见表1;岩体选择HJC本构方程描述,SPH方法中各参数见文献[11];混凝土材料选择C&K材料描述[18],钢筋和锚杆选用率相关塑性随动强化模型(*MAT_PLASTIC_KINEMATIC),参数见文献[19]。TNT炸药采用*MAT_HIGH_EXPLOSIVE_BURN材料模型结合*EOS_JWL状态方程描述[20]:

表1 钨合金弹壳的材料参数Table 1 Material parameters of tungsten alloy cartridge case

(3)

式中:P为爆轰压力;V为相对体积;E0表示初始体积内能;A、B、R1、R2和ω为炸药常数。炸药参数取值见表2。

表2 炸药的材料参数Table 2 Material parameters of explosives

钢筋与混凝土、锚杆与岩石通过关键字*CONSTRAINED_LAGRANGE_IN_SOLID实现相互协调变形、共同受力。

3 模型验证

由于缺乏同等当量试验实测资料,采用侵彻深度、爆炸峰值压力计算公式验证,以确保模型建立和参数确定的正确性。文献[21]中通过对比多组国内外岩体侵彻数据发现,Bernard公式(C)对于不同RQD值的岩体均具有较高计算精度。Bernard公式(C)的无量纲形式表示为:

(4)

图3为弹头侵彻深度时程曲线和速度时程曲线,由图可见:t=15 ms时起爆炸药;弹头在侵入岩体时速度值降为0。模拟中弹壳壁薄,刚度比实际刚度小,装药强度小,另外从弹头开始装药,与实际装药有差异,装药强度比弹体强度低,在侵彻过程中,弹体在钻地后会发生挤压变形,比实际的刚性弹侵彻下的阻力变大,因此模拟侵彻深度与理论计算值相比较小,误差为1.9%,模拟精度较高。

图3 弹头速度、深度时程曲线Fig.3 Time history curve of warhead velocity and depth

计算爆炸远区岩体峰值压力判断模型有效性其压力峰值σp可表示为[22]:

(5)

式中:R为爆炸传播距离;c为介质中冲击波波速;f为耦合系数;ρ为岩石密度;W为装药质量;n为介质衰减系数。

选择爆炸远区岩体单元作为研究对象,式(5)计算的单元峰值压力值与模拟值对比见表3。由表3可知,模拟单元压力值均大于计算峰值,因为侵彻过程中已经形成初始应力,最大误差为7.9%。

表3 数值模拟与计算公式结果对比Table 3 Comparison of results of numerical simulation and empirical formula

以上验证说明,本构模型选择和参数确定合理,数值模拟结果可信。

4 结果与分析

4.1 侵彻爆炸联合作用过程分析

Mises强度理论认为材料内一点的Mises等效应力达到某材料极限强度时破坏,其取值为:

(6)

式中:σ1,σ2,σ3分别指第一、第二、第三主应力。

图4为计算模型的Mises应力云图,由于采用cm-g-ms单位制计算,图中所示单位为105MPa。由图4可见,在侵彻爆炸过程中,不断有岩体粒子向上向外抛出,侵爆产生的高压气团急速压缩周围岩体,并在岩体中形成空腔。

图4 侵彻爆炸不同时刻Mises应力云图Fig.4 Mises stress nephogram at different moments of penetration and explosion

图4(a)初始t=8 ms纯侵彻时,弹体应力较大,岩体中仅形成细条状弹坑,此时弹头处最大等效应力为1.436 GPa,超过屈服强度1.2 GPa,弹头与炸药衔接处沿径向张开,这是因为炸药与弹壳材料强度数值上相差较大,强冲击作用下两种材料界面处出现应力集中,包裹炸药处合金截面较小使得应力较大而最先屈服,并逐渐发生变形;图4(b)为t=15 ms时炸药爆炸,弹尾出现较大应力变化,细条状弹坑伴随弹腹扩张向两侧扩展并呈现纺锤状,这与岩体内柱形装药爆炸初期爆腔形状一致,弹头继续向深部侵彻并发生较大变形,这是因为受爆炸压力和弹体反作用力叠加作用,大部分弹体单元等效应力超过屈服强度;图4(c)为t=29.5 ms时,钻地弹主体基本解体,零星破片随大量岩体向外抛出,应力波继续向深部岩体传播,纺锤状应力云图扩大,原因是近地表岩体由于受到侵彻初始挤压影响,密度变大,应力波传播较快,在横向传播距离更长,出现上宽下窄的类纺锤形;图4(d)为t=46 ms时,抛掷物变少,被覆结构直墙拱底应力发生变化,表明前驱应力波已穿过地下结构。SPH粒子与六面体网格界面处应力云图连续,进一步证明了岩体损伤近区SPH算法与FEM算法耦合有效,参数设置合理。

4.2 封堵作用对侵彻爆炸的影响

通常FEM算法模拟侵彻爆炸过程是分段进行的,即首先模拟弹体侵彻岩体过程,侵彻结束后再在侵蚀算法形成弹坑处安置等当量炸药,模拟爆炸过程,该方法存在以下问题:一是对于侵彻效应的模拟采用侵蚀算法,即单元应力或应变达到破坏准则即被删除,如图5所示,网格的删除造成了质量不守恒,整个过程只能大致符合侵彻爆炸的物理过程;二是由于侵彻通道单元的删除,形成了侵彻孔,对爆炸应力波有卸载作用,类似于卸压爆破,忽略了弹壳及侵彻中坍落岩体对弹坑的封堵作用,分段模拟的方法形成的Mises应力云图如图6所示,形成等效峰值应力为246.8 MPa。实际上,坚硬弹体的高速冲击对岩体进行了快速而剧烈的挤压,弹头近区岩石密度升高,对应单元即使挤压破碎失效也同样会对弹体产生握裹阻力。SPH-FEM耦合算法模拟侵彻爆炸过程是连续的,可以考虑侵彻后弹体尾部失效散落岩体对爆炸应力波的封堵作用和弹头近区岩石的握裹力。考虑封堵作用下爆心等效峰值应力为437 MPa,比卸荷爆破增加了77%。对比图4(c)和图6(a)可以发现,侵彻爆炸连续作用下爆腔形状较卸压爆破模拟的椭圆爆腔形状反映了弹头形状的挤压变形和近地端半径更大,更符合弹体侵彻爆炸的实际情况。

4.3 直墙拱破坏模式分析

拱体混凝土材料采用K&C模型,定义比例损伤因子δ来表示其损伤情况,其表示为:

(7)

式中:λ为材料损伤函数值;λm为材料初始损伤量。LS-Prepost中用标记为“plastic strain”的变量云图代替比例损伤因子表示的损伤云图,如图7所示。

图7(a)t=7.5 ms初始侵彻阶段时,由于拱体与围岩无空隙,拱顶外侧直接承受上部压力荷载,损伤明显;图7(b)t=13 ms时,损伤由拱顶向下蔓延,拱顶损伤穿透结构,拱肩损伤由外向内延伸,破坏区域占1/2,尚未形成贯穿破坏;图7(c)t=35 ms时,爆炸产生应力波到达结构,损伤进一步加剧,损伤传递至直墙底部,由于底板与结构拱分离浇筑,垂直振动时,结构产生竖向变形,应力波产生的结构振动对底板影响不大,未发生明显破坏;图7(d)t=47 ms时,底板与直墙相邻处发生损伤,中心无损伤变化;图7(e)t=50 ms时通过结构底部局部放大图可以看到,底板与围岩间形成了一定空隙。

为进一步了解底板与围岩间产生空隙原因,选取图7(e)单元为研究对象,分别绘制直墙上、中和下部侧方对应围岩单元R1、R2和R3的X向位移时程曲线,如图8(a);底板1/2、1/4和边缘处下部对应围岩单元R4、R5和R6的Y向位移时程曲线,如图8(b);下述X、Y轴正方向如图7(e)。由图8(a)可知,R1、R2和R3单元都有先向X轴正向移动至峰值后逐渐减小,随后产生负向位移,且随着单元Y负向深度的增加,X轴向位移峰值逐渐变小且时间延后,这表明在应力波作用下,受成拱效应和振动效应影响,围岩先受压向X轴正向位移,而后由于围岩反向振动发生回弹导致其负向位移;由图8(b)可知,底板下的围岩先向Y轴负向移动至峰值后减小,随后又向Y轴正向移动,且距离直墙脚越近位移越大,这说明上部荷载主要由直墙脚处向下传递,较大的压力使得脚部发生较大位移,由于底板与主体拱结构分离浇筑导致底板下围岩受动载相对较小,底板中心几乎没有位移。综上所述,围岩与拱结构的相互作用使得底板受到向内挤压从而发生翘曲,底板与直墙采用分离式浇筑,塑性铰弯矩变为挤压力,可以减轻地板的损伤,直墙脚处双向受压损伤较大。

图8 围岩位移时程曲线Fig.8 Time history curve of surrounding rock displacement

为从单元应力变化情况具体了解上述损伤破坏传递情况,图9给出了拱顶、拱肩、拱脚和直墙底部单元Von Mises等效应力时程曲线。

图9 结构关键点等效应力时程曲线Fig.9 Time history curve of effective stress of key points of structure

在结构破坏过程中,由于侵彻爆炸联合作用应力波的叠加,各部位应力出现明显双峰或多峰情况。拱顶在爆炸主导阶段出现3个峰值,分别为9.01 MPa、13.49 MPa和16.05 MPa,拱肩应力峰值较小,均不大于7.71 MPa,但在破坏过程中多次出现;拱脚在爆炸主导阶段出现最大峰值(18.36 MPa);直墙底部应力峰值出现在破坏末期,峰值较小(6.90 MPa)。整个过程荷载作用情况复杂,在炸药起爆前(t=15.0 ms)各部位应力逐渐增加,但最大值不大于峰值的50%,说明纯侵彻作用对结构毁伤效果有限;但后期侵彻爆炸两者应力波的复杂叠加,使结构内部不仅有高应力产生,而且峰值多次出现,对结构抗性提出更高要求。