含初始缺陷的复合固体推进剂力学性能

封 涛, 许进升, 韩 龙, 陈 雄, 周长省

复合固体推进剂作为一种含能复合材料被广泛用作运载火箭及各种战略战术导弹的动力能源，其力学性能很大程度上影响着导弹的生存能力及作战能力。早期对复合固体推进剂的研究大多基于连续介质力学，通过实验从唯象的角度得到推进剂宏观力学响应的本构关系[1-4]，但这种方法无法对推进剂内在结构变化机理进行有效研究。近年来随着计算机性能的提高，国内外众多学者开始着手从细观角度出发，建立了复合固体推进剂细观颗粒的填充模型，对推进剂的力学性能进行了定量的分析，并取得了一定的成果。

Matous等[5]最先通过自主开发的Rocpack软件生成了固体推进剂代表性体积单元，并在颗粒与基体之间的界面层设置了粘接单元来模拟界面脱粘损伤的产生与发展。常武军等[6]通过实验证实了颗粒/基体界面是推进剂的薄弱环节，其脱粘是造成宏观应力-应变曲线非线性的重要原因之一，获得界面确切的力学性能参数对推进剂的细观研究工作至关重要。张兴高等[7]采用扫描电镜，并结合界面化学原理，准确测量了AP及HTPB黏合剂基体的接触角和表面能参数，得到了界面性能与老化之间的有效关系。Zhi等[8-10]在研究推进剂的细观损伤和非线性力学性能时，用双线性内聚力模型近似表征推进剂颗粒/基体界面的脱粘损伤，同时还研究了颗粒的体积分数、粒径大小、位置的随机分布对推进剂细观损伤及宏观力学性能的影响。韩龙等[11]研究了复合固体推进剂细观界面性能随拉伸速率的变化规律，通过参数反演的分析算法得到了优化后的界面参数值，该值能准确预测不同加载速率下推进剂的宏观力学行为。

以往的研究中多认为推进剂为黏合剂基体、增强颗粒和界面相组成的三相复合材料，认为其内部细观结构完好，没有考虑到真实情况下推进剂生产过程中存在的初始缺陷的影响。事实上，复合固体推进剂细观组成与混凝土极为相似，二者经常被用来作比较研究，受混凝土研究中的启发[12-13]，本研究认为，由于工艺限制，推进剂细观颗粒/基体界面处存在着包裹不完全的初始缺陷，并对含缺陷的界面力学性能进行了定义，研究当界面存在缺陷时推进剂的宏观力学行为，得到推进剂的初始模量及拉伸强度随界面缺陷含量的变化规律。

1 细观建模及参数获取

复合固体推进剂是由分散相颗粒、中间相界面和连续相基体组成的高填充比三相复合材料。要从细观层面上研究其力学性能，首要的是建立一个能符合推进剂真实细观结构的颗粒填充模型，然后赋予模型相应的力学参数和边界条件，模拟推进剂在载荷作用下的力学响应。

1.1 细观建模

1.1.1 代表性体积单元生成

颗粒增强复合材料细观建模主要有基于蒙特卡罗算法[14]的随机模拟方法和基于分子动力学算法[15]的碰撞模拟方法，由于复合固体推进剂的颗粒填充体积分数通常较大，可以达到70%～80%，再综合考虑到算法实现的效率、成熟度以及有限元前处理的工作量等相关因素，本研究选用分子动力学方法建立复合固体推进剂的颗粒填充模型。

表1为复合固体推进剂的典型配方。图1为推进剂典型配方中颗粒粒径分布。

表1 复合固体推进剂典型配方Table 1 Typical formula of composite solid propellant

结合表1所示推进剂配方及图1所示的颗粒分布规律，应用分子动力学方法对复合固体推进剂颗粒填充模型建模，如图2所示。由于颗粒的填充体积分数较高，颗粒间大小差异很大，在这种模型上直接生成网格进行有限元计算是非常困难的。为了简化计算模型，参考文献[16]通过Mori-Tanaka等效的方法将Al颗粒对推进剂力学性能的影响等效到复合基体中，在计算时只单独考虑AP颗粒的影响。

1.1.2 边界条件

基于周期性假设认为，推进剂在外载荷的作用下，其内部的应力场和应变场存在连续性和周期性，因此在利用有限元方法预测复合材料固有属性参数时，需要对代表性体积单元（representative volume element，RVE）施加周期性边界条件以满足宏观上的连续均质假设，并确保得到合理的细观应力、应变分布。周期性边界可以表示为：

式中：L为RVE单元的尺寸；和为施加在边界上的位移载荷。

由于周期边界条件须在周期性网格的基础上方能实现，在实际应用中较为不便，为了能在有限元软件ABAQUS中既准确又简便地模拟推进剂的单轴拉伸力学行为，本研究选用均匀位移边界条件，如图2（a）所示。均匀位移边界条件下，RVE在受载变形过程中各边始终保持平直，图2（a）中虚线为变形后的RVE边界位置。上边界所给的均匀位移载荷为0.1667 mm/s，对应的应变率为0.003333 s–1，同时给右边界施加一个限制x方向位移一致的耦合约束。文献[17]研究表明，使用均匀位移边界与周期性边界计算的结果不存在明显的差异，并且随着RVE尺寸的增加误差渐渐趋近于零。

1.1.3 网格划分

本研究是基于推进剂二维颗粒填充模型展开，并视为平面应变问题进行数值求解。在有限元计算中，AP颗粒被视为弹性体，其模量相对基体材料大得多，近似认为在受载时不变形，于是在网格划分时对其采用四节点的平面应变单元CPE4，并在界面附近处将网格细化以提高计算的收敛性。

在有限元分析中，除平面应力问题外，需采用杂交单元来模拟不可压缩材料（泊松比为0.5）或近似不可压缩材料（泊松比大于0.475）的响应。HTPB基体材料为近似不可压缩材料，在受载荷时体积保持恒定，不随应力应变状态而发生改变，在计算时无法利用节点处的位移插值函数来获得单元的压应力，而杂交单元毋须通过节点位移插值，通过自身设定的额外自由度便可直接获取单元的压应力，其偏应力及偏应变的获取则仍然依赖于节点的位移场计算。因此对基体材料的网格划分选用平面应变4节点四边形线性积分的杂交单元CPE4H，模型网格划分如图2（b）所示。

在颗粒与基体界面处设置粘接单元（cohesive element）模拟颗粒与基体的粘接特性，并沿界面的径向设置扫掠网格，为保证粘接单元的收敛，在计算时设置较小的初始载荷步并对粘接单元设定一定的黏性。通过在界面处定义摩擦罚函数来施加通用接触算法以防止界面脱粘后颗粒与基体可能会发生重叠或相互渗透。

1.2 组分参数获取

复合推进剂对应的组分参数为黏合剂基体的松弛模量和固体填充颗粒的弹性模量及泊松比。其中颗粒的相关参数依据文献[5]选取，AP：弹性模量32447 MPa，泊松比 0.1433；Al：弹性模量 68300 MPa，泊松比0.33。

对推进剂基体材料的等速率拉伸应力响应，本研究从线黏弹性理论出发，结合对基体材料松弛实验[18]测得的应力响应，将松弛模量以Prony级数形式进行拟合，具体拟合表达式为[19]：

式中：为平衡模量；和分别为第i个Maxwell单元的模量和松弛时间；t为时间。松弛曲线拟合结果见图3。由此可得HTPB基体的黏弹性松弛行为描述，进而得到基体基于松弛模量的应力-应变本构关系[1]：

式中：为基体应力响应，为应变率。

2 界面缺陷单元定义

内聚力模型（cohesive zone model，CZM）建立在弹塑性断裂力学的基础上，最早由Dugdale[20]和Barenblatt[21]在研究脆性材料的断裂时提出来的，其基本思想是将界面理想化为具有一定粘接强度等力学特性的无厚度面，并通过牵引力-位移法则（tractionseparation law，TSL）定义粘接单元的力学响应，从而表征整个界面处的损伤起始和演化形式。

2.1 双线性内聚力模型

在有限元ABAQUS中，内聚区采用一层厚度近似为0的内聚力单元表示，将内聚力单元嵌入到传统单元之间使其上下表面与相邻单元连接，外载荷引起的材料损伤只发生在内聚力单元中，周围单元不受影响。通过定义内聚力单元的TSL本构关系，以实现内聚力单元在外力作用下损伤的起始与演化。双线型损伤内聚力模型由于结构简单，在较复杂的工程环境下也能很好地表征内聚区损伤演化等力学行为，被广泛地应用于各种研究之中，典型的双线型损伤内聚力模型[22]示意图如图4所示。

其满足的牵引力-位移法则如下：

式中：下标，t 分别代表界面的法向和切向；为内聚应力；为界面的分离位移；和分别表示法向和切向的内聚强度，且；是临界位移，表示内聚应力达到内聚强度时对应界面的分离位移；为定义的损伤变量。

二维界面的内聚力单元中存在两个应力分量和位移分量，载荷作用初期，界面应力和位移为线性关系，当应力达到界面损伤初始准则后，界面处产生损伤，其应力和位移关系不再线性变化，本研究采用式（6）中的二次应力作为损伤起始准则：

损伤变量的定义如式（6）：

式中：为失效位移，表示界面分离位移的最大值，此后界面破坏失效，失去承载能力。

采用式（7）断裂能准则判断界面是否破坏失效：

式中：，分别为界面的法向和切向的内聚能，为界面参数（取值为2）。

2.2 初始缺陷定义

复合固体推进剂是典型的多相混合非均匀材料，在其制备过程中由于材料特性及制备工艺限制，使得推进剂中会随机分布着大量的初始缺陷，包括：在机械混合时无法严格控制真空环境，使得有少许空气进入黏合剂基体而导致孔洞和气泡；在加入固体填充颗粒后搅拌至混合均匀过程中，颗粒间因相互碰撞导致的颗粒破碎；更多的缺陷是在固化降温过程中，由于颗粒和基体导热性能差异使得二者变形不一致而导致的颗粒/基体界面缺陷。这些初始缺陷的存在降低了推进剂的结构强度和刚度，严重降低了其力学性能。因此需要对含初始缺陷的复合固体推进剂颗粒填充模型进行数值计算，从细观尺度探究初始缺陷的存在对推进剂力学性能的影响。

以往有关推进剂的宏观力学实验证实[6]，颗粒/基体界面作为推进剂结构中最薄弱的一环，在载荷作用下界面的脱粘是导致推进剂失效破坏的根本原因。本研究只考虑推进剂颗粒/基体界面处的初始缺陷，并注意到推进剂在固化降温过程中形成的应力集中与颗粒间距有关这一事实，对界面缺陷作如下假设：

（1）对于某一给定位置的颗粒，在其周边两倍最大颗粒粒径的同心圆域内，计算其他颗粒与该颗粒之间的距离，如果该间距小于设定的最小距离，即时认为在该颗粒界面处存在初始缺陷，本研究中，的取值为最小填充颗粒的粒径。

（2）初始缺陷在两颗粒外公切线区域内的界面上均匀分布，如图5所示，采用弱化的界面单元来表征初始缺陷，定义缺陷系数为，则含缺陷的界面力学性能为，，。

（3）设界面初始缺陷含量为，其定义为，其中为单个界面缺陷单元周长；为所有界面单元总周长。

2.3 界面力学参数确定

细观界面力学性能参数的取值直接影响到数值计算结果的准确性，目前还不能通过实验技术测得其确切值，对于参数的选取更多的是依据文献和经验来预估，但是这样得到的模型参数往往并不精确，因此，众多学者提出了基于Hooke-Jeeves的参数反演方法[23]，其基本原理是通过不断地调整模型参数的取值来改变仿真结果，最终使仿真结果与实验值的误差小于规定值的过程。

为反演需要，对复合固体推进剂进行单轴拉伸实验，不考虑推进剂力学性能的温度相关性和应变率相关性，实验在室温下进行，采用的拉伸速率为10 mm/min，对应的应变率为 0.003333 s–1，得到推进剂的工程应力-应变曲线，结果如图6所示。

界面力学模型中主要有三个模型参数，界面临界位移、内聚强度和破坏位移。经研究发现，改变临界位移只会影响仿真结果中应力-应变曲线的线性上升段的斜率，曲线的下降段不会受到影响。内聚强度的取值主要影响应力-应变曲线的峰值载荷区域，究其原因，内聚强度是弹性段与损伤段的分界点，当应力达到内聚强度之后，损伤起始，界面单元力学性能开始发生衰减。破坏位移表示的是界面单元失效时界面层上下端的分离位移，其对应力-应变曲线下降段有显著的影响，取值越大曲线下降段上载荷值越大。针对这三种参数对应力-应变曲线的不同影响特性，确定了分步反演的思路，构建表征实验曲线与仿真曲线重合度目标函数：

当目标函数的值小于选定容差时，认为当前界面参数仿真所得结果与实验结果吻合，界面参数的初始值选取及反演优化后的值列于表2中。

表2 界面内聚参数Table 2 Cohesive parameters of interface

2.4 界面参数验证

在2.3节已通过推进剂的单轴拉伸实验反演出粘结完好的颗粒/基体界面的内聚参数，但还需要进一步验证参数的准确性。考虑到固体火箭发动机装药实际工作过程中承载的复杂性及多样性，采用多阶段的加载方式进行验证实验。首先将推进剂试样以10 mm/min的速率拉伸至10%应变，接着松弛1 min，再以10 mm/min的速率对试样加载至20%应变。同时，运用本研究构建的细观模型和相关界面参数，对该实验过程进行有限元仿真分析，输出验证曲线，验证实验曲线与仿真曲线，如图7所示。由图7可以看出，实验结果与仿真结果吻合度较高，曲线的大体走势基本一致，虽然两者仍存在一定的偏差，但在一定程度上仿真结果能够反映真实实验过程中的载荷变化情况，说明本研究通过反演优化算法获得的界面参数是准确可靠的，能用以描述推进剂颗粒与基体界面的粘接力学性能。

3 算例分析

采用商业有限元软件ABAQUS，建立了4种颗粒含量均为0.72具有不同分布的复合固体推进剂细观颗粒填充模型（图8），并对所建立的模型进行数值计算。为满足对比分析的需要，模型（a）暂未考虑初始缺陷的影响，模型（b），（c），（d）均基于2.2节中的假设定义初始缺陷，并计算得到其对应的缺陷含量分别为10%，20%和30%。限于篇幅，这里仅取界面缺陷含量为20%的模型（c）为例，图9给出了相应初始缺陷分布及载荷作用下应力-应变云图。

从图9中可以看出，当整体应变为5%时，由于推进剂组分材料属性的不同以及颗粒之间的相互作用，导致推进剂内部应力和应变分布很不均匀。颗粒作为弹性体，其模量相对于基体和界面大得多，在拉伸过程中，颗粒的变形较小，基体的变形较大，这就使得颗粒/基体界面成为整个推进剂中的薄弱环节，在载荷作用下很容易就会产生界面脱粘现象。尤其当颗粒/基体界面处存在着初始缺陷时，即使很小的外力也能使这些弱界面处产生很明显的脱粘形貌，从而产生微孔洞。

当应变达到15%，颗粒/基体界面处的变形进一步增加，在定义的缺陷界面处开始出现脱粘，形成明显的微孔洞。究其原因，推进剂中颗粒与基体形变的差异以及颗粒的相互作用很大程度决定了推进剂内部的应力和应变分布，进而影响了颗粒与基体粘接界面处的脱粘过程。由于在定义的缺陷界面处其局部应变会先达到相应的失效位移，导致界面破坏失去承载能力，形成明显的孔洞。在以往文献中多提出大颗粒附近界面会较先产生脱粘，其根本原因是在界面性能相同的前提下，较大颗粒界面处的局部应变更大，最先达到界面的失效位移致使界面失效，这与本研究得出的弱界面最先产生脱粘失效虽看似不同，但两者的本质是一样的，都是局部应力超过了当前的界面强度引起的。

整体应变持续增加达到50%，这时推进剂中因界面脱粘产生的微孔洞处就开始出现明显的裂纹直至破坏。因为随着载荷的进一步增加，界面进一步扩展，颗粒与基体的分离程度也不断增加，不仅定义的弱界面几乎都已经发生脱粘失效，在较大颗粒处粘接完好的界面也因局部应变达到失效位移而失效。界面完全脱粘后，颗粒不再承受应力，从而导致在因脱粘形成的裂隙尖端附近存在明显应力集中现象，慢慢地这些由于界面脱粘而导致基体应力集中处的微孔洞逐渐汇合，最终使基体撕裂导致裂纹产生，进而造成推进剂的宏观破坏。由图9中真实推进剂拉断的电镜扫描图可见，红框标记的脱粘形貌与仿真结果吻合较好，验证了本研究方法的可靠性。

将仿真结果中模型上边界各节点在加载方向上的约束反力求和即为上边界所受的拉力，经处理可得推进剂仿真所得的工程应力-应变曲线，将其与前节的实验曲线的平均值进行对比，如图10所示。

比较图9中推进剂拉伸应力云图和图10应力-应变曲线可以看出：在较小应变（）时，推进剂处于线弹性段，应力随着应变的增加表现出线性增长，在该阶段推进剂内部颗粒/基体界面处虽有较大应变，但仍未达到界面的失效位移，即使存在缺陷的颗粒界面处也无脱粘现象。继续加载，当应变达到15%时，由图9可知，此时推进剂内部缺陷界面处存在着密集的脱粘，界面承载能力下降，在宏观应力-应变曲线上此时开始出现非线性，进入“应力平台区”。载荷继续增加到图9所示的50%应变时，这时即便当初粘接完好的颗粒基体界面也产生了严重的脱粘，脱粘处的颗粒丧失承载能力，黏合剂基体完全承担外载荷，而后推进剂基体被不断拉长，最终因微孔洞的汇合造成基体撕裂，在很短时间内推进剂宏观破坏。

比较图10中不含缺陷的模型仿真曲线与推进剂的宏观实验曲线，可以看出：经过反演优化获得的仿真曲线与实际实验曲线吻合度较高，但同时也存在着一定的差距。这可能是因为：（1）本研究的推进剂细观仿真是建立在二维平面应变基础上的，这与真实推进剂三维结构应力状态存在着一定的偏差；（2）所建立的细观模型中所有固体颗粒形貌均是规则的圆形，实际推进剂中可能含有一些不规则的颗粒，相比于规整的球体这些颗粒附近更容易出现应力集中，在载荷作用下更容易产生脱粘现象，这也是为什么实验应力-应变曲线始终在仿真曲线下方的原因。

比较图10中不同界面缺陷含量的推进剂仿真结果曲线，不难发现：（1）界面缺陷的存在使得推进剂的初始模量降低了，并且缺陷含量越高其模量下降得越剧烈，这是因为含有缺陷的界面传递载荷的能力相对较低，在相同的小应变下，其对应的应力也较低，缺陷含量越高则表现的更明显；（2）界面缺陷的存在降低了推进剂的拉伸强度，这主要是因为存在缺陷的颗粒/基体界面在拉伸过程中会更早出现脱粘，造成推进剂应力下降，缺陷含量愈高则下降愈多。图11为初始模量、拉伸强度与缺陷含量关系。初步探索表3中推进剂的初始模量及拉伸强度随缺陷含量的变化规律，根据图11中散点分布，认为其可能为线性或指数关系，在尝试选用线性拟合时发现当缺陷含量为1时，推进剂的初始模量和拉伸强度均为负，这违背了物理规律，因此本研究选用指数函数拟合其关系，相应的拟合结果见图11。

表3 推进剂力学性能Table 3 Mechanical properties of composite solid propellant

4 结论

（1）构建了复合固体推进剂的细观颗粒填充模型，并基于合理的假设定义了含缺陷界面处的力学性能参数，通过仿真计算发现缺陷界面处会更早达到损伤条件发生脱粘，缺陷的存在加速了推进剂的断裂失效的进程。

（2）比较了不同缺陷含量对应的推进剂应力-应变曲线，发现界面缺陷的存在降低了推进剂的力学性能，缺陷含量越高则力学性能下降越明显，推进剂的初始模量及拉伸强度随缺陷含量的增加呈指数下降的趋势。

［1］JUNG G D, YOUN S K. A nonlinear viscoelastic constitutive model of solid propellant[J]. International Journal of Solids & Structures, 1999, 36（25）: 3755-3777.

［2］XU F, ARAVAS N, SOFRONIS P. Constitutive modeling of solid propellant materials with evolving microstructural damage[J]. Journal of the Mechanics & Physics of Solids, 2008, 56（5）: 2050-2073.

［3］ 王哲君, 强洪夫, 王广, 等. 固体推进剂力学性能和本构模型的研究进展[J]. 含能材料, 2016, 24（4）: 403-416.

（WANG Z J, QIANG H F, WANG G, et al. Review on the mechanical properties and constitutive models of solid propellants[J]. Chinese Journal of Energetic Materials,2016, 24（4）: 403-416.）

［4］BECKER E B. Constitutive modeling of high-elongation solid propellants[J]. Journal of Engineering Materials and Technology, 1992, 114: 111.

［5］MATOUŠ K, INGLIS H M, GU X, et al. Multiscale modeling of solid propellants: from particle packing to failure[J]. Composites Science & Technology, 2007, 67（7/8）: 1694-1708.

［6］ 常武军, 鞠玉涛, 王蓬勃. HTPB推进剂脱湿与力学性能的相关性研究[J]. 兵工学报, 2012, 32（3）: 261-266.

（CHANG W J, JU Y T, WANG P B. Research on correlation between dewetting and mechanical property of HTPB propellant[J]. Acta Armamentarii, 2012, 32（3）:261-266.）

［7］ 张兴高, 张炜, 芦伟, 等. HTPB推进剂填料/基体界面粘结性能老化特性研究[J]. 含能材料, 2009, 17（3）: 269-273.

（ZHANG X G, ZHANG W, LU W, et al. Aging characteristics of interfacial adhesive property of filler/binder matrix for HTPB propellant[J]. Chinese Journal of Energetic Materials, 2009, 17（3）: 269-273.）

［8］ZHI S J, SUN B, ZHANG J W. Multiscale modeling of heterogeneous propellants from particle packing to grain failure using a surface-based cohesive approach[J]. Acta Mechanica Sinica, 2012, 28（3）: 746-759.

［9］ 职世君, 曹付齐, 申志彬, 等. 复合固体推进剂颗粒脱湿损伤参数反演[J]. 推进技术, 2016, 37（10）: 1977-1983.

（ZHI S J, CAO F Q, SHEN Z B, et al. Parameters inversion of particle dewetting damage of composite solid propellants[J]. Journal of Propulsion Technology, 2016,37（10）: 1977-1983.）

［10］职世君, 孙冰, 张建伟. 基于表面粘结损伤的复合固体推进剂细观损伤数值模拟[J]. 推进技术, 2013, 34（2）:273-279.

（ZHI S J, SUN B, ZHANG J W. Numerical simulation of solid sropellant mesoscopic damage using surfacebased cohesive approach[J]. Journal of Propulsion Technology, 2013, 34（2）: 273-279.）

［11］韩龙, 许进升, 周长省. HTPB/IPDI复合固体推进剂细观界面率相关参数的反演识别研究[J]. 含能材料, 2016,24（10）: 928-935.

（HAN L, XU J S, ZHOU C S. Inverse identification of the rate-dependent micro interface parameters of HTPB/IPDI composite propellant[J]. Chinese Journal of Energetic Materials, 2016, 24（10）: 928-935.）

［12］王娟, 管巧艳, 孔宇田, 等. 界面初始缺陷对混凝土单轴抗拉强度的影响[J]. 水利水电科技进展, 2016, 36（5）:46-49.

（WANG J, GUAN Q Y, KONG Y T, et al. Influences of initial interfacial defects on uniaxial tensile strength of concrete[J]. Advances in Science and Technology of Water Resources, 2016, 36（5）: 46-49.）

［13］杜修力, 揭鹏力, 金浏. 考虑初始缺陷影响的混凝土梁动态弯拉破坏模式分析[J]. 工程力学, 2015（2）: 74-81.

（DU X L, JIE P L, JIN L. Dynamic flexural-tensile failure mode analysis of concrete beam with initial defect[J].Engineering Mechanics, 2015（2）: 74-81.）

［14］FISH J. Multiscale modeling and simulation of composite materials and structures[J]. Multiscale Methods in Computational Mechanics, 2011, 55: 215-231.

［15］HUMPHREYS D D, FRIESNER R A, BERNE B J. A multiple-time-step molecular dynamics algorithm for macromolecules[J]. Journal of Physical Chemistry, 2002,98（27）: 6885-6892.

［16］马昌兵. 复合固体推进剂细观结构建模及其力学行为数值模拟[D]. 西安: 第二炮兵工程学院, 2011.

（MA C B. Research on mesostructured modeling and mechanical behaviors numerical simulation of composite solid propellant[D]. Xi’an: The Second Artillery Engineering University, 2011.）

［17］PIERARD O. Micromechanics of inclusion-reinforced composites in elasto-plasticity and elasto-viscoplasticity:modeling and computation[D]. Leuven: University Catholique de Louvain, 2006.

［18］常武军. 复合固体推进剂细观损伤及其数值仿真研究[D]. 南京: 南京理工大学, 2013.

（CHANG W J. Research on microstructural damage and its numerical simulation method for composite solid propellant[D]. Nanjing: Nanjing University of Science and Technology, 2013.）

［19］XU J, HAN L, ZHENG J, et al. Finite element implementation of a thermo-damage-viscoelastic constitutive model for hydroxyl-terminated polybutadiene composite propellant[J]. Mechanics of Time-Dependent Materials,2017: 1-19.

［20］DUGDALE D S. Yielding of steel sheets containing slits[J]. Journal of the Mechanics & Physics of Solids,1960, 8（2）: 100-104.

［21］BARENBLATT G I. The mathematical theory of equilibrium cracks in brittle fracture[J]. Advances in Applied Mechanics, 1962, 7: 55-129.

［22］ELICES M, GUINEA G V, GÓMEZ J, et al. The cohesive zone model: advantages, limitations and challenges[J]. Engineering Fracture Mechanics, 2002, 69（2）: 137-163.

［23］JEEVES T A. “Direct Search” solution of numerical and statistical problems[J]. Journal of the Acm, 1961, 8（2）:212-229.