中国指挥与控制学会会刊 
现代战争中,破障是研究大规模登陆作战无法回避的重要课题,在破障行动开始前,登陆作战一方需要对通路开辟的类型、数量、位置和长宽等进行精确的计算,并提出合理的破障需求[1-3]。阻绝墙是反登陆作战的重要手段,可以阻碍进攻方前进路线,迟滞对方行动。目前,破障弹是破除阻绝墙、开辟通路的主要武器装备。研究者从多个角度对破障弹进行了研究,比如王刚进行了聚能装药载荷下混凝土破坏行为的实验研究[4],郑平泰进行了射流侵彻混凝土介质数值模拟及影响因素研究[5],刘燕燕进行了线性聚能射流作用下岩石裂纹扩展数值模拟[6]等。以上研究对混凝土与岩石板材料的毁伤效果进行了实验与仿真研究,但是缺少对特定作战条件下的不同材料、不同厚度阻绝墙的毁伤效能具体数值的研究。
本文拟在破障弹的实战运用环境下,通过分析阻绝墙破障弹的毁伤机理,建立阻绝墙破障弹毁伤效能的数值计算模型,分别对混凝土靶板与岩石靶板进行破障弹毁伤效能仿真评估,为提升作战运用效能提供新的技术方法。
1 试验材料与数据分析方法
1.1 阻绝墙破障弹基本参数
阻绝墙破障弹是基本参数为130 mm×1 000 mm的锥形弹,卵形头部尺寸约为200 mm,总质量为17.6 kg,载药当量约为10 kg,外观如图1所示。炸药爆炸后可在阻绝墙中形成冲击波或者应力波,使介质产生裂缝,随后爆生气体产生,继而膨胀,其膨胀过程延长了应力波的作用时间,且爆生气体挤入裂缝产生气刃作用,致使裂缝扩展,最终完成破障任务。
1.2 阻隔墙类型及特点
根据材料的不同,常用的阻绝墙可以分为C30混凝土、钢质结构、岩石结构等,其弹性模量、泊松比、屈服强度、本构关系等存在较大差别,因而其爆轰破坏模式、毁伤效果也有较大区别。C30混凝土与常用岩石板的材料力学性能参数见表1。
表1 阻绝墙材料参数Tab.1 Material parameters of blocking wall |
| 材料 | 密度/ (kg/m3) | 弹性模量/ MPa | 屈服强度/ MPa |
|---|---|---|---|
| C30混凝土 | 2 360 | 28 000 | 40 |
| 岩石 | 2 500 | 40 000 | 100 |
2 破障弹毁伤效能仿真数值计算模型
传统的有限元法是(隐式)结构有限元,而毁伤效应分析常用的有限元法是显式动力有限元,两者的基本原理和求解过程大体相同[7],主要区别是结构有限元需要由单元刚度矩阵集成总体刚度矩阵,而显式动力有限元不形成总体刚度矩阵,避免了矩阵求逆,减小了计算量。此外,显式动力有限元法采用中心差分时间积分,这是一种条件稳定算法,采用的时间步长必须小于某个临界值,否则算法不稳定[8]。有限元法将求解区域划分为有限个在节点处相连接的小单元,单元可以具有不同的形状,以不同的连接方式进行组合,很容易划分复杂的几何形状。常用的有限元单元有四边形单元、六面体单元和四面体单元等[9]。下面以常用的八节点六面体单元为例介绍有限元法的空间离散[10]。单元内任意点的坐标可以用节点坐标插值表示为
$x_{i}(\xi, \eta, \zeta, t)=\sum_{j=1}^{8} \varphi_{j}(\xi, \eta, \zeta) x_{i}^{j}(t), i=1,2,3$
式中:(ξ,η,ζ)为单元的自然坐标; (t)为t时刻第j节点的坐标值。插值函数(形函数)ϕi (ξ,η,ζ)为
$\varphi_{j}=\frac{1}{8}\left(1+\xi \xi_{j}\right)\left(1+\eta \eta_{j}\right)\left(1+\zeta \zeta_{j}\right), j=1,2,3$
式中:(ξj,ηj,ζj)为单元第j节点的自然坐标。式(1)用矩阵可表示为
$\boldsymbol{x}(\xi, \eta, \zeta, t)=\boldsymbol{N} \boldsymbol{x}^{e}=\left[x_{1}, x_{2}, x_{3}\right]$
其中,差值矩阵为
$\boldsymbol{N}(\xi, \eta, \zeta)=\left[\begin{array}{l}\varphi_{1}, 0,0, \varphi_{2}, 0, \cdots, 0,0 \\0, \varphi_{1}, 0,0, \varphi_{2}, \cdots, \varphi_{8}, 0 \\0,0, \varphi_{1}, 0,0, \cdots, 0, \varphi_{8}\end{array}\right]_{3 \times 24}$
设整个求解域离散为n个单元,虚功方程近似为
$\delta \Pi=\sum_{m=1}^{n} \delta \Pi_{m}=0$
写成矩阵形式,可以得到
$\begin{array}{c}\sum_{m=1}^{n} \delta x^{e \mathrm{~T}}\left(\int_{v_{m}} \rho N^{\mathrm{T}} N \mathrm{~d} V \ddot{x}^{e}+\int_{v_{m}} B^{\mathrm{T}} \sigma \mathrm{~d} V-\right. \\\left.\int_{v_{m}} \rho N^{\mathrm{T}} f \mathrm{~d} V-\int_{S_{2}^{m}} N^{\mathrm{T}} t \mathrm{~d} S\right)=0\end{array}$
式中,B为6×24阶的应变矩阵。Cauchy应力矢量为
$\boldsymbol{\sigma}_{T}=\left[\boldsymbol{\sigma}_{x}, \boldsymbol{\sigma}_{y}, \boldsymbol{\sigma}_{z}, \boldsymbol{\sigma}_{x y}, \boldsymbol{\sigma}_{y z}, \boldsymbol{\sigma}_{z x}\right]$
由动量守恒方程得
$\boldsymbol{M} \ddot{\boldsymbol{x}}(t)=-\boldsymbol{F}(x, \dot{x})+\boldsymbol{P}(x, t)$
式中:M为总体质量矩阵; 为总体节点加速度矢量;P为总体节点荷载矢量,由节点集中力、面力、体力等形成;F由单元应力场的等效节点力矢量(或称应力散度)组集而成:
$\boldsymbol{F}=\sum_{m=1}^{n} \int_{v_{m}} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{\sigma} \mathrm{~d} V$
结构大变形时伴随存在大转动,在本构方程中与应变率对应的应力率必须是关于刚体转动具有不变性的客观张量,在连续介质力学中,Jaumann应力率σΔ就是这种张量。Jaumann应力率 的表达式为
$\boldsymbol{\sigma}_{i j}^{\Delta}=\dot{\boldsymbol{\sigma}}_{i j}-\boldsymbol{\sigma}_{i k} \boldsymbol{\Omega}_{j k}-\boldsymbol{\sigma}_{j k} \boldsymbol{\Omega}_{j k}$
式中:Ωij为旋率张量,即
$\boldsymbol{\Omega}_{i j}=\frac{1}{2}\left(\frac{\partial \dot{x}_{i}}{\partial x_{j}}-\frac{\partial \dot{x}_{j}}{\partial \dot{x}_{i}}\right)$
显式有限元算法一般采用单点高斯(Gauss)积分,可以极大地节省数据存储量和运算次数,但单点高斯积分可能引起零能模式,单点高斯积分时单元变形的沙漏模态被丢失,即它对单元应变能的计算没有影响,故称零能模式。在动力响应计算时,沙漏模态将不受控制,导致计算结果出现数值震荡。显式有限元算法一般采用沙漏黏性阻尼来控制零能模式。考虑体积黏性阻尼力的影响,运动方程可以改写为
$\boldsymbol M \ddot{\boldsymbol{x}}(t)=\boldsymbol{P}-\boldsymbol{F}+\boldsymbol{H}-\boldsymbol{C} \dot{\boldsymbol{x}}$
采用中心差分法求时间积分,即
$\ddot{\boldsymbol{x}}\left(t_{n}\right)=\boldsymbol{M}^{-1}\left[\boldsymbol{P}\left(t_{n}\right)-\boldsymbol{F}\left(t_{n}\right)+\boldsymbol{H}\left(t_{n}\right)-\boldsymbol{C}\left(t_{n-1 / 2}\right)\right]$
$\dot{\boldsymbol{x}}\left(t_{n+1 / 2}\right)=\dot{\boldsymbol{x}}\left(t_{n-1 / 2}\right)+\frac{1}{2}\left(\Delta t_{n-1}+\Delta t_{n}\right) \ddot{\boldsymbol{x}}\left(t_{n}\right)$
$\boldsymbol{x}\left(t_{n+1}\right)=\boldsymbol{x}\left(t_{n}\right)+\Delta t_{n} \dot{\boldsymbol{x}}\left(t_{n+1 / 2}\right)$
式中:tn-1/2=(tn+tn-1)/2,tn+1/2=(tn+tn+1)/2,Δtn-1=(tn-tn-1),Δtn=(tn+1-tn)。 (tn)、 (tn+1/2)、x(tn+1)分别为tn时刻的节点加速度矢量、tn+1/2时刻的节点速度矢量和tn+1时刻的节点坐标矢量。
由于采用集中质量矩阵,上述运动方程的求解是非耦合的,不需要集成总体矩阵,节省了存储空间和计算时间,但是显式中心差分法是有条件稳定的。
3 破障弹毁伤效能仿真
3.1 破障弹毁伤效能仿真参数设定
在阻绝墙毁伤仿真分析中,一定当量炸药爆炸效应下,不同材料的阻绝墙力学响应差异较大,且由于化学反应过程极为复杂,研究人员难以使用数学模型对其毁伤效果进行描述,因此,根据实际情况,本文选定的阻绝墙材料为混凝土、岩石靶板。人们建立了相同尺寸的阻绝墙有限元模型,墙体高度约为10×1.5 m,为评估其毁伤效果,厚度分别取为50 mm、100 mm、200 mm以及500 mm。爆炸中心(即着靶点)设在墙体正面中心处,炸药模型采用JWL模型,当量设置为10 kg。起爆时间定义为瞬时起爆,不同计算单元中DAMAGE参量达到1或者单元自动删除则视为有效毁伤,据此进行毁伤范围计算。
3.2 破障弹毁伤效能仿真结果分析
图2为单发破障弹对混凝土材质阻绝墙的毁伤效果示意,在实际使用中,每发破障弹的毁伤孔洞半径至少达到625 mm,才能算作爆破成功。计算结果表明,不同厚度混凝土板毁伤方式基本一致,主要差异在于毁伤程度。当厚度为50 mm、100 mm以及200 mm时,弹体爆炸可沿中心产生损伤孔洞,其尺寸分别约为829 mm、737 mm以及497 mm,见图2(a,b,c)。厚度增加,空洞半径则降低,这是因为厚度越大,阻绝墙可以吸收更多爆炸能量,其抵抗变形及失效的能力也越强。图2(d)所示,在厚度达到500 mm时,阻绝墙不产生明显的孔洞,即破障弹难以直接炸穿此厚度及以上的混凝土阻绝墙。除此之外,图中红色部分代表该区域损伤程度达到失效标准,对应真实情况则为该部分阻绝墙体产生明显的破坏失稳,与其余基体的连接强度大幅降低,使用铲车及其他清障设备可以更轻松地进行拆除,4种工况塑性失稳范围大致相同,约为1 800 mm。结果表明,在进行低于200 mm厚度的混凝土阻绝墙毁伤时,1发破障弹即可满足预期效果,破损墙体可利用工程机械装备进行清理,即可开拓宽度为2 m的作战通道。当厚度达到500 mm时,单发破障弹作用下阻绝墙受破坏较大,但墙体材料并未完全脱落,因此,可补射1发破障弹并进行障碍清理,即可开拓相应的作战通道。
4 结束语
本文评估了单发破障弹对不同材质、不同厚度阻绝墙的毁伤效果,结果表明,相比于岩石靶板,破障弹对混凝土砌体有较好的毁伤效果。在进行低于200 mm厚度的混凝土阻绝墙毁伤时,单发破障弹即可满足预期效果,当厚度达到500 mm时,需要补射破障弹并进行障碍清理。根据大量实弹射击基础上提出的别列赞公式计算,本文研究的破障弹对混凝土的有效毁伤厚度约303.2 mm,对岩石靶板的有效毁伤厚度约373.1 mm,仿真结果与试验结果基本一致。当岩石结构厚度达到500 mm时,需要补射多发破障弹才能打通作战通路。