1 引言

传统的炸药配方设计以试验为主，需要大量的炸药性能试验确定配方，设计成本高，周期长。近年来，虚拟设计技术发展迅速，其本质是以计算机支持的模拟仿真技术为前提，在产品设计阶段实时预测产品的性能，从而提高产品的一次成功率。虚拟设计技术已在炸药配方设计领域得到了应用，包括爆轰性能、力学性能等炸药基本性能的预估^[1-7]。

在炸药有效弹性性能的虚拟仿真设计方面，国内主要采用分子动力学(MD)方法^[2]，而国外则较多地采用细观尺度数值仿真的方法，Biswajit Banerjee和E.M.Mas采用细观尺度仿真模型计算了PBX炸药的有效弹性性能^[4-6]，S.Ravi Annapragada计算了B炸药的弹性模量和泊松比等有效弹性性能^[7]。国外的研究中，细观尺度计算模型的建立多基于实际炸药粒径的分布，即仿真预估前已经得到了炸药样品，并进行了粒度测试。然而在炸药配方虚拟设计阶段没有炸药试样，无炸药粒径分布的准确数据，不能采用上述细观尺度仿真的方法来预估配方的有效弹性性能。

鉴于此，我们建立了B炸药有效弹性性能的细观尺度仿真模型，并根据文献[7]的数据对其进行了验证。在此基础上，进一步研究了粒径分布对计算结果的影响。结果表明：即使没有准确的炸药粒度分布数据，也能准确预估B炸药的有效弹性性能，可以应用于炸药配方虚拟设计。

2 仿真模型

数值仿真采用有限元方法，通过构建细观尺度的炸药代表体积元(RVE)，计算静力加载下炸药RVE的弹性响应，进而求出弹性模量和泊松比。

2.1 RVE建模方法^[5-7]

在平面应变假设下建立二维的计算模型，RVE为正方形，炸药分散相颗粒的形状作圆形简化，圆形颗粒填充在炸药连续相基体中。RVE建模分为以下步骤：

(1) 设定RVE的边长L，生成RVE的基体模型。L的确定主要根据炸药分散相颗粒的最大粒径d_max，并结合计算硬件条件进行设定。根据相关研究，L可以小于d_max，也可以为d_max的几倍。L越大，包含的分散相颗粒越多，计算规模越大，因此L不宜太大，在满足计算精度的前提下尽量小一些。

(2) 在得到分散相颗粒的粒径分布的前提下，按照粒径从大到小的顺序在基体上填充颗粒。首先填充直径最大的颗粒，填充过程中满足以下条件：(a)圆心坐标随机生成；(b)颗粒之间不能重叠；(c)通过设置各颗粒之间的距离的上下限来控制颗粒的团聚程度；(d)颗粒不能全部位于基体以外，如果颗粒越过了基体边界，那么须满足循环边界条件；(e)预先设定颗粒的体积分数，等达到预先设定的体积分数时，完成第一批颗粒的填充。

(3) 按照上述填充方法，依次进行不同粒径颗粒的填充，当所有的颗粒的体积分数等于炸药分散相体积分数时，填充过程结束。

按照上述思路，在ANSYS程序中采用APDL参数化语言编制了相应的建模程序，用于生成RVE几何模型。

RVE应具有循环对称边界，设已经填充的颗粒的圆心坐标为(x_i，y_i)，半径为r_i，待填入颗粒的圆心坐标为(x_j，y_j)，半径为r_j，颗粒间距为Δl。图 1为不越过边界的情形，对此进行颗粒填充时须满足：

$ \begin{array}{l} \sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l,\\ \left( {i = 1,2 \cdots \cdots ,j = 2,3, \cdots \cdots i + 1} \right) \end{array} $

(1)

图 2为越过右侧边界的情形，对此要在模型左侧补充另一个颗粒以保证循环边界条件，此时须满足：

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {\sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l}\\ {\sqrt {{{\left( {{x_i} - a - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l} \end{array}} \right.\\ \left( {i = 1,2, \cdots \cdots ,j = 2,3, \cdots \cdots i + 1} \right) \end{array} $

(2)

颗粒越过单侧边界的其他情况可以参照(2)式处理。如果颗粒同时越过两个边界(如图 3所示)，此时应满足：

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {\sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l}\\ {\sqrt {{{\left( {{x_i} - a - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l}\\ {\sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} + a - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l}\\ {\sqrt {{{\left( {{x_i} - a - {x_j}} \right)}^2} + {{\left( {{y_i} + a - {y_j}} \right)}^2}} > {r_i} + {r_j} + \Delta l} \end{array}} \right.\\ \left( {i = 1,2, \cdots \cdots ,j = 2,3, \cdots \cdots i + 1} \right) \end{array} $

(3)

同样地，颗粒同时越过两侧边界的其他情况可以参照(3)式处理，在此不再赘述。

2.2 边界条件和载荷

计算过程包括设置材料参数、边界条件和载荷，划分网格和执行计算。计算所需的材料参数是混合炸药各单质组分的有效弹性性能参数。

根据2.1所述方法建立典型RVE模型(见图 4)，其边长为0.75 mm，RVE中颗粒共有三种，直径分别为0.25，0.13，0.05 mm，体积分数分别为0.276, 0.183, 0.141，全部颗粒的总体积分数为0.6。图 4同时给出了计算时RVE边界条件和载荷，并满足下式：

$ \left\{ {\begin{array}{*{20}{l}} {{u_x}\left( {x,y} \right) = 0\;\;on\;\;x = 0}\\ {{u_y}\left( {x,y} \right) = 0\;\;on\;\;y = 0}\\ {\varepsilon = 0.01\;\;\;\;\;\;\;\;on\;\;y = L} \end{array}} \right. $

(4)

式中, u_x和u_y分别为x方向和y方向的位移。

2.3 有效弹性性能计算公式

按照下式计算炸药的弹性模量和泊松比^[5]:

$ \left\{ \begin{array}{l} E_{{\rm{eff}}}^{2{\rm{D}}} = {E_x} = \frac{{{F_x}}}{{{u_x}}}, v_{{\rm{eff}}}^{2{\rm{D}}} =-\frac{{{u_y}}}{{{u_x}}}\\ {u_x} = \frac{1}{m}\sum\limits_{i = 1}^m {{u_{xi}}}, {u_y} = \frac{1}{n}\sum\limits_{i = 1}^m {{u_{yi}}} \end{array} \right. $

(5)

式中，E_eff^2D和ν_eff^2D分别为二维弹性模量和泊松比，F_X为x方向反应力，m，n分别是模型在x和y方向的节点数，u_xi和u_yi分别是边界节点在x方向和y方向的位移。

由于上述结果是基于平面应变假设而得到，所以还要将二维计算结果换算到三维当中，即：

$ \left\{ \begin{array}{l} v_{{\rm{eff}}}^{3{\rm{D}}} = \frac{{v_{{\rm{eff}}}^{2{\rm{D}}}}}{{1 + v_{{\rm{eff}}}^{2{\rm{D}}}}}\\ E_{{\rm{eff}}}^{3{\rm{D}}} = E_{{\rm{eff}}}^{2{\rm{D}}}\left( {1-{{\left( {v_{{\rm{eff}}}^{2{\rm{D}}}} \right)}^2}} \right) \end{array} \right. $

(6)

式中，E_eff^3D和ν_eff^3D为三维弹性模量和泊松比。

3 仿真模型验证

选取B炸药验证上述仿真模型的精度。如同文献[7]，本研究中选取三种粒径RDX颗粒：50, 130, 250 μm，三种RDX颗粒在B炸药中所占的体积分数分别为13.3%, 17.2%, 26.0%。本研究中RVE的边长为1 mm，B炸药二维RVE模型如图 5所示。划分网格时单元的边长为5 μm，共生成90152个计算单元。计算时输入的材料参数见表 1。

图 6是RVE模型的等效应力分布云图，根据图 6可以看出，RDX颗粒的等效应力比TNT基体的等效应力大，这是由于RDX的弹性模量较大，导致在相同应变下的应力较大。根据计算得到边界处的节点位移和反应力，基于公式(5)和(6)，计算得到了B炸药的弹性模量和泊松比，并与文献[7]中给出的试验结果进行了对比，结果见表 2。由表 2可见，弹性模量和泊松比的计算值与实验值的误差分别为8.26%和3.44%，都在10%以内，满足工程预估精度要求。这说明采用RVE方法可以准确计算出混合炸药的有效弹性性能，只要能够获得准确的混合炸药单质组分的力学性能，就能通过RVE方法准确预测混合炸药的有效弹性性能。

4 粒径分布对有效弹性性能的影响

本研究RVE的建模过程中使用了炸药颗粒的粒径分布数据，然而在炸药配方虚拟设计阶段，尚无法得到炸药的粒径分布。因此，需要研究颗粒的粒径分布对于计算结果的影响。仍以B炸药作为研究对象，建模时不考虑RDX颗粒的详细粒径分布，假设各个RDX颗粒的粒径都相等。设定RVE的边长为1 mm，RVE中RDX颗粒的数量分别为8, 12, 20, 40个，分别得到了相应的RVE模型，如图 7。

根据公式(5)和(6)计算出了B炸药的弹性模量和泊松比，结果见表 3。由表 3可知：在RVE建模时不考虑RDX的粒径分布，同样能得到准确的有效弹性性能结果，计算值与试验值的误差不超过10%。这说明基于有限元的RVE方法并不受限于炸药颗粒的粒径分布，也就是说，在混合炸药配方设计阶段，即使没有掌握详细的粒径分布数据，也可以采用RVE方法预估混合炸的有效弹性性能。

5 结论

采用ANSYS程序APDL参数化编程语言编制了炸药细观尺度RVE建模程序，并以B炸药为研究对象计算了其有效弹性性能，计算结果与文献结果吻合良好。在此基础上进一步研究了炸药颗粒粒径分布对于计算结果的影响，结果表明：不详细考虑炸药颗粒的粒径分布，预估结果也有较高的精度，有效弹性性能的计算误差不超过10%，这说明本文的方法可以应用到炸药配方仿真设计阶段的力学性能预估。