论文部分内容阅读
摘要： 为研究复合固体推进剂松弛模量的衰减， 依据复合固体推进剂组分建立单胞体模型。 根据复合固体推进剂5%定应变下的松弛模量， 利用有限元计算方法和均匀化理论对基体材料的松弛模量进行参数反演。 在此基础上， 研究了复合固体推进剂单胞体模型在不同定应变水平下的松弛模量， 并计算了单胞体模型的拉伸过程。 根据Boltzmann叠加原理及广义胡克定律计算了单胞体模型拉伸过程中松弛模量的变化。 结果表明， 复合固体推进剂松弛模量有明显的应变相关性， 并且在拉伸过程中模量的衰减与载荷历程相关。
关键词： 复合固体推进剂； 松弛模量； 细观力学； 数值模拟
中图分类号： V435文献标识码： A文章编号： 1673-5048（2015）02-0050-05
Abstract： To research the decrement of solid propellant relaxation modulus， the single cell model is established based on the components of a solid propellant. According to the solid propellant relaxation modulus which obtained under 5% strain load， the binder relaxation modulus is inversed by the FEM and homogenization theory. On this basis， the relaxation modulus of the solid propellant are studied when the solid propellant under different strain loads， and the drawing process of the solid propellant is computed. The variation of relaxation modulus in the drawing process of the solid propellant is obtained based on Boltzmann superposition principle and the generalized Hooke law. The results show that relaxation modulus of the composite solid propellant is dependent on its strain loads， and the decrement of relaxation modulus is dependent on the loading history in the drawing process.
Key words： composite solid propellants； relaxation modulus； mesomechanics； numerical simulation
0引言
随着世界军事科技的发展， 导弹性能的不断提高， 对固体火箭发动机性能的要求越来越高。 如何在保证固体火箭发动机装药结构完整性的前提下，
尽可能地提高发动机的装药量， 一直是固体火箭发动机设计过程中需要解决的关键问题。 要解决该问题， 首先必须要展开固体推进剂力学性能的研究。 复合固体推进剂是一种多相（基体、 增强相、 界面相等）复合材料， 其力学性能受增强相的体积分数及其组分材料性质的影响比较大。 因此， 基于细观力学理论研究复合固体推进剂非线性力学性能， 不仅可以直观地反映出组分材料的影响， 还可以对其宏观非线性本构关系的研究提供理论支撑。
复合固体推进剂是一种高填充比的复合材料。 近些年来， 随着计算机性能的大幅度提高， 计算细观力学得到了迅速发展。 Matous等[1-2]在复合固体推进剂颗粒与基体之间的界面层设置了粘结单元（Cohesive element）模拟了固体颗粒和基体之间损伤的产生及发展。 Tan等[3-4]通过数字图象等相关技术获得了高能炸药PBX9501紧凑拉伸试样裂尖周围的应力场及位移场， 利用扩展的MoriTanaka方法对试验结果做了均匀化处理， 并结合试验研究得到的颗粒与基体之间的非线性粘结模型， 对含不同尺寸颗粒的细观模型进行了数值计算。 国内一些学者研究了复合固体推进剂的颗粒夹杂模型的建模方法[5-7]。 在此基础上， 有的学者采用有限元法对复合固体推进剂进行直接数值模拟， 对推进剂内部界面脱粘过程进行了有限元分析[8]。 有的学者根据数值仿真结果， 结合细观力学方法， 如： MoriTanaka方法或改进的MoriTanaka方法， 研究了固体推进剂的模量、 界面脱粘对固体推进剂力学性能的影响等[9-11]， 或采用多步法， 通过将基体与部分颗粒均质化为一种混合物， 计算出较复杂的固体推进剂的有效模量[12]。
本文在之前研究[13-14]的基础上， 通过建立复合固体推进剂单胞体模型， 根据复合固体推进剂松弛模量试验值对基体材料的松弛模量进行参数反演， 结合数值计算方法对复合固体推进剂松弛模量的衰减特性进行了研究。
1物理模型和计算方法
1.1 颗粒夹杂模型及组分材料属性
复合固体推进剂是一种典型的颗粒增强复合材料， 由于颗粒粒径的分布特征及颗粒的随机分布， 使得其微结构特征非常复杂， 如图1所示。
从图1中可以看出， 复合固体推进剂颗粒堆积模型可以更真实地反映其微结构特征， 但这种模型建模过程复杂， 并且由于颗粒之间距离较近， 使得网格密度较大， 计算量较大， 这在后续工作中将逐步展开研究。 本文旨在研究复合固体推进剂松弛模量应变的衰减特性， 因此首先从单胞体模型出发进行研究。 根据复合固体推进剂的固体含量确定AP颗粒的体积分数， 本文所研究的复合固体推进剂AP颗粒体积分数为65%， 单胞体模型如图2所示。 由之前研究[13]可知， 复合固体推进剂的松弛特性主要来源于基体材料， 在AP颗粒的增强作用下， 复合固体推进剂的松弛特性与基体材料基本一致， 区别主要体现在瞬时模量E0的增大。 因此， 基于图2所示模型， 根据复合固体推进剂的松弛模量对基体材料的松弛模量进行参数反演。 首先确定基体瞬时模量E0的上、 下限值E10和E20， 并根据已确定的上、 下限值进行有限元计算， 结合细观力学方法， 得出基体瞬时模量E0的上、 下限值所对应的复合固体推进剂瞬时模量E1P0和E2P0； 然后对比复合固体推进剂的实际瞬时模量EP0， 根据EP0在E1P0和E2P0之间的线性分布关系， 在E01和E20之间确定下一步计算时新的基体瞬时模量E01′， 取E10和E20的平均值为新的E02′， 再次进行有限元计算； 反复进行迭代， 最后确定基体的松弛模量。 假设AP颗粒为弹性体， 取其弹性模量和泊松比分别为： E=32 450 MPa， ν=0.143 3[1]。
1.2均匀化方法
当复合固体推进剂单胞体模型受载荷作用时， 其体积平均应力和平均应变可由式（1）～（2）计算：
当AP颗粒体积分数为65%时， 图2模型中颗粒之间距离很小， 基体材料部分若采用六面体网格， 在颗粒之间临近部分， 网格易产生畸变。 因此， 对于AP颗粒可采用六面体网格， 基体材料部分采用四面体网格。 根据有限元高斯积分原理， 分别计算单个六面体单元和四面体单元的平均应力、 应变。 模型的体积平均应力和应变可以通过单个六面体网格和四面体网格的平均应力和平均应变来确定， 即分别表示第m个四面体单元的平均应力和平均应变； nHexij和nHexij分别表示第n个六面体单元的平均应力和平均应变； NTet和NHex分别表示四面体单元和六面体单元的数量； VTetm， VHexn和VRVE分别表示第m个四面体单元， 第n个六面体单元及整个单胞体的体积。
2计算结果及讨论
复合固体推进剂松弛模量一般是根据GJB-770B在定应变5%下测量的。 根据复合固体推进剂定应变为5%时的松弛模量， 结合均匀化方法对基体材料的松弛模量进行参数反演， 获得基体材料的松弛模量。 定应变为5%时单胞体模型拉伸方向上的应力分布如图3所示。 从图中可以看出， 应力较大的区域主要集中在颗粒间距小的区域， 随着时间的增大， 应力有明显的松弛效应。
为研究复合固体推进剂松弛模量的应变相关性， 假设基体的松弛模量不变， 应变小于8%时忽略界面损伤的影响， 分别计算当定应变为0.01%， 1%， 2%， 3%， 4%， 5%， 6%， 7%， 8%时复合固体推进剂的松弛模量， 如图4所示。 从图中可以明显看出， 当定应变较小时所计算的松弛模量较大， 随着定应变值的增大， 复合固体推进剂松弛模量逐渐减小。 所得结论与züpek[15]根据不同定应变（0.5%， 1%， 2%， 5%， 10%和15%）下测得的松弛模量试验结果得到的结论一致。
分别对单胞体模型及复合固体推进剂均质材料模型的拉伸过程进行数值计算。 其中， 复合固体推进剂均质材料模型采用不同定应变下的松弛模量。 拉伸时间分别取0.05 s和100 s， 最终的真实应变均为7.16%。 计算结果如图5和图6所示。 从图中可以看出， 采用不同定应变水平下的松弛模量的拉伸曲线与单胞体模型的拉伸曲线有明显差异， 只有当真实应变约小于3%时， 所选取定应变值越低松弛模量的计算结果与单胞体模型的计算结果越接近。
因此， 可根据单胞体模型的拉伸曲线， 选取一定的时间增量Δt， 求解出t时刻时的松弛模量E（t）， 对比初始时刻的松弛模量曲线， 从而获得瞬时模量E0的变化。 为验证该方法及确定合适的时间增量Δt， 以复合固体推进剂定应变为5%时的松弛模量为输入条件， 拉伸时间为0.05 s， 取Δt分别为0.001 s， 0.000 5 s和0.000 1 s计算复合固体推进剂的松弛模量， 并与输入的松弛模量进行对比， 相对误差如图7所示。 从图中可以看出， 相对误差随时间的增大而增大， 这是因为随着计算时间的增加， 在t～t+Δt时刻内， 时间t之前累计载荷的松弛效应逐渐增大， 而在计算的过程中忽略了Δt时刻内松弛效应的影响。 虽然相对误差随时间的增大而增大， 但数值仍然较小， 说明了该计算方法可行。 从不同时间增量Δt的计算结果来看， 当Δt从0.001 s缩小10倍后， 相对误差产生了震荡， 并且减小Δt时相对误差降低的幅度不大。 输入条件不变， 拉伸时间改为100 s， Δt取0.001 s， 计算的松弛模量与理论值对比如图8所示。 当拉伸时间为100 s时， 相对误差为0.83%。 采用该方法可以根据复合固体推进剂的拉伸曲线计算松弛模量的变化。
3总结
（1） 由于复合固体推进剂材料的颗粒夹杂效应， 通过仿真计算可知其松弛模量具有明显的应变相关性， 即在不同的定应变水平下测得的试验结果会存在差异， 随着定应变水平的增大， 测得的松弛模量逐渐减小。 计算结果所得结论与züpek[15]试验结果一致。
（2） 根据单胞体模型不同拉伸速率下， 不考虑界面损伤时的拉伸曲线， 计算了拉伸过程中复合固体推进剂松弛模量的变化。 计算结果表明， 在拉伸过程中， 复合固体推进剂的松弛模量有明显的衰减现象， 并且衰减程度与载荷历程相关。
参考文献：
[1] Matous K， Inglis H M， Gu X F， et al. Multiscale Damage Modeling of Solid Propellants： Theory and Computational Framework[C]//41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit， AIAA 2005-4347， 2005： 1-14. [2] Matous K， Geubelle P H. Multiscale Modeling of Particle Debonding in Reinforced Elastomers Subjected to Finite Deformation[J]. International Journal for Numerical Methods in Engineering， 2006， 65（2）： 190-223.
[3] Tan H， Liu C， Huang Y， et al. The Cohesive Law for the Particle/Matrix Interfaces in High Explosives[J]. Journal of the Mechanics and Physics of Solids， 2005， 53（8）： 1892-1917.
[4] Tan H， Huang Y， Liu C， et al. The Uniaxial Tension of Particulate Composite Materials with Nonlinear Interface Debonding[J]. International Journal of Solids and Structures. 2007， 44（6）： 1809-1822.
[5] 史佩， 李高春， 李昊. 复合固体推进剂细观力学模型研究[J]. 计算机仿真， 2007， 24（5）： 21-24.
[6] 赵玖玲， 田先斌. 颗粒复合材料代表性体元并行建模算法研究[J]. 计算机仿真， 2010， 27（1）： 46-49.
[7] 刘著卿， 李高春， 邢耀国. 复合固体推进剂细观损伤扫描电镜实验及数值模拟[J]. 推进技术， 2011， 32（3）： 412-416.
[8] 曲凯， 张旭东， 李高春. 基于内聚力界面脱粘的复合固体推进剂力学性能研究[J]. 火炸药学报， 2008， 31（6）： 77-81.
[9] 李高春， 邢耀国， 王玉峰. 基于细观力学的复合固体推进剂模量预估方法[J]. 推进技术， 2007， 28（4）： 441-444.
[10] 李高春， 邢耀国， 戢治洪， 等. 复合固体推进剂细观界面脱粘有限元分析[J]. 复合材料学报， 2011， 28（3）： 229-235.
[11] 刘承武， 阳建红， 陈飞. 改进的Mori-Tanaka法在复合推进剂非线界面脱粘中的应用[J]. 固体火箭技术， 2011， 34（1）： 67-70.
[12] 马昌兵， 强洪夫， 武文明， 等. 颗粒增强复合材料有效弹性模量预测的多步法[J]. 固体力学学报， 2010， 31（S1）： 12-16.
[13] Zhi Shijun， Sun Bing， Zhang Jianwei. 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.
[14] 职世君， 孙冰， 张建伟. 基于表面粘结损伤的复合固体推进剂细观损伤数值模拟[J]. 推进技术， 2013， 34（2）： 273-279.
[15] Sebnem züpek. Constitutive Equations for Solid Propellants [D]. University of Texas at Austin， 1997.
关键词： 复合固体推进剂； 松弛模量； 细观力学； 数值模拟
中图分类号： V435文献标识码： A文章编号： 1673-5048（2015）02-0050-05
Abstract： To research the decrement of solid propellant relaxation modulus， the single cell model is established based on the components of a solid propellant. According to the solid propellant relaxation modulus which obtained under 5% strain load， the binder relaxation modulus is inversed by the FEM and homogenization theory. On this basis， the relaxation modulus of the solid propellant are studied when the solid propellant under different strain loads， and the drawing process of the solid propellant is computed. The variation of relaxation modulus in the drawing process of the solid propellant is obtained based on Boltzmann superposition principle and the generalized Hooke law. The results show that relaxation modulus of the composite solid propellant is dependent on its strain loads， and the decrement of relaxation modulus is dependent on the loading history in the drawing process.
Key words： composite solid propellants； relaxation modulus； mesomechanics； numerical simulation
0引言
随着世界军事科技的发展， 导弹性能的不断提高， 对固体火箭发动机性能的要求越来越高。 如何在保证固体火箭发动机装药结构完整性的前提下，
尽可能地提高发动机的装药量， 一直是固体火箭发动机设计过程中需要解决的关键问题。 要解决该问题， 首先必须要展开固体推进剂力学性能的研究。 复合固体推进剂是一种多相（基体、 增强相、 界面相等）复合材料， 其力学性能受增强相的体积分数及其组分材料性质的影响比较大。 因此， 基于细观力学理论研究复合固体推进剂非线性力学性能， 不仅可以直观地反映出组分材料的影响， 还可以对其宏观非线性本构关系的研究提供理论支撑。
复合固体推进剂是一种高填充比的复合材料。 近些年来， 随着计算机性能的大幅度提高， 计算细观力学得到了迅速发展。 Matous等[1-2]在复合固体推进剂颗粒与基体之间的界面层设置了粘结单元（Cohesive element）模拟了固体颗粒和基体之间损伤的产生及发展。 Tan等[3-4]通过数字图象等相关技术获得了高能炸药PBX9501紧凑拉伸试样裂尖周围的应力场及位移场， 利用扩展的MoriTanaka方法对试验结果做了均匀化处理， 并结合试验研究得到的颗粒与基体之间的非线性粘结模型， 对含不同尺寸颗粒的细观模型进行了数值计算。 国内一些学者研究了复合固体推进剂的颗粒夹杂模型的建模方法[5-7]。 在此基础上， 有的学者采用有限元法对复合固体推进剂进行直接数值模拟， 对推进剂内部界面脱粘过程进行了有限元分析[8]。 有的学者根据数值仿真结果， 结合细观力学方法， 如： MoriTanaka方法或改进的MoriTanaka方法， 研究了固体推进剂的模量、 界面脱粘对固体推进剂力学性能的影响等[9-11]， 或采用多步法， 通过将基体与部分颗粒均质化为一种混合物， 计算出较复杂的固体推进剂的有效模量[12]。
本文在之前研究[13-14]的基础上， 通过建立复合固体推进剂单胞体模型， 根据复合固体推进剂松弛模量试验值对基体材料的松弛模量进行参数反演， 结合数值计算方法对复合固体推进剂松弛模量的衰减特性进行了研究。
1物理模型和计算方法
1.1 颗粒夹杂模型及组分材料属性
复合固体推进剂是一种典型的颗粒增强复合材料， 由于颗粒粒径的分布特征及颗粒的随机分布， 使得其微结构特征非常复杂， 如图1所示。
从图1中可以看出， 复合固体推进剂颗粒堆积模型可以更真实地反映其微结构特征， 但这种模型建模过程复杂， 并且由于颗粒之间距离较近， 使得网格密度较大， 计算量较大， 这在后续工作中将逐步展开研究。 本文旨在研究复合固体推进剂松弛模量应变的衰减特性， 因此首先从单胞体模型出发进行研究。 根据复合固体推进剂的固体含量确定AP颗粒的体积分数， 本文所研究的复合固体推进剂AP颗粒体积分数为65%， 单胞体模型如图2所示。 由之前研究[13]可知， 复合固体推进剂的松弛特性主要来源于基体材料， 在AP颗粒的增强作用下， 复合固体推进剂的松弛特性与基体材料基本一致， 区别主要体现在瞬时模量E0的增大。 因此， 基于图2所示模型， 根据复合固体推进剂的松弛模量对基体材料的松弛模量进行参数反演。 首先确定基体瞬时模量E0的上、 下限值E10和E20， 并根据已确定的上、 下限值进行有限元计算， 结合细观力学方法， 得出基体瞬时模量E0的上、 下限值所对应的复合固体推进剂瞬时模量E1P0和E2P0； 然后对比复合固体推进剂的实际瞬时模量EP0， 根据EP0在E1P0和E2P0之间的线性分布关系， 在E01和E20之间确定下一步计算时新的基体瞬时模量E01′， 取E10和E20的平均值为新的E02′， 再次进行有限元计算； 反复进行迭代， 最后确定基体的松弛模量。 假设AP颗粒为弹性体， 取其弹性模量和泊松比分别为： E=32 450 MPa， ν=0.143 3[1]。
1.2均匀化方法
当复合固体推进剂单胞体模型受载荷作用时， 其体积平均应力和平均应变可由式（1）～（2）计算：
当AP颗粒体积分数为65%时， 图2模型中颗粒之间距离很小， 基体材料部分若采用六面体网格， 在颗粒之间临近部分， 网格易产生畸变。 因此， 对于AP颗粒可采用六面体网格， 基体材料部分采用四面体网格。 根据有限元高斯积分原理， 分别计算单个六面体单元和四面体单元的平均应力、 应变。 模型的体积平均应力和应变可以通过单个六面体网格和四面体网格的平均应力和平均应变来确定， 即分别表示第m个四面体单元的平均应力和平均应变； nHexij和nHexij分别表示第n个六面体单元的平均应力和平均应变； NTet和NHex分别表示四面体单元和六面体单元的数量； VTetm， VHexn和VRVE分别表示第m个四面体单元， 第n个六面体单元及整个单胞体的体积。
2计算结果及讨论
复合固体推进剂松弛模量一般是根据GJB-770B在定应变5%下测量的。 根据复合固体推进剂定应变为5%时的松弛模量， 结合均匀化方法对基体材料的松弛模量进行参数反演， 获得基体材料的松弛模量。 定应变为5%时单胞体模型拉伸方向上的应力分布如图3所示。 从图中可以看出， 应力较大的区域主要集中在颗粒间距小的区域， 随着时间的增大， 应力有明显的松弛效应。
为研究复合固体推进剂松弛模量的应变相关性， 假设基体的松弛模量不变， 应变小于8%时忽略界面损伤的影响， 分别计算当定应变为0.01%， 1%， 2%， 3%， 4%， 5%， 6%， 7%， 8%时复合固体推进剂的松弛模量， 如图4所示。 从图中可以明显看出， 当定应变较小时所计算的松弛模量较大， 随着定应变值的增大， 复合固体推进剂松弛模量逐渐减小。 所得结论与züpek[15]根据不同定应变（0.5%， 1%， 2%， 5%， 10%和15%）下测得的松弛模量试验结果得到的结论一致。
分别对单胞体模型及复合固体推进剂均质材料模型的拉伸过程进行数值计算。 其中， 复合固体推进剂均质材料模型采用不同定应变下的松弛模量。 拉伸时间分别取0.05 s和100 s， 最终的真实应变均为7.16%。 计算结果如图5和图6所示。 从图中可以看出， 采用不同定应变水平下的松弛模量的拉伸曲线与单胞体模型的拉伸曲线有明显差异， 只有当真实应变约小于3%时， 所选取定应变值越低松弛模量的计算结果与单胞体模型的计算结果越接近。
因此， 可根据单胞体模型的拉伸曲线， 选取一定的时间增量Δt， 求解出t时刻时的松弛模量E（t）， 对比初始时刻的松弛模量曲线， 从而获得瞬时模量E0的变化。 为验证该方法及确定合适的时间增量Δt， 以复合固体推进剂定应变为5%时的松弛模量为输入条件， 拉伸时间为0.05 s， 取Δt分别为0.001 s， 0.000 5 s和0.000 1 s计算复合固体推进剂的松弛模量， 并与输入的松弛模量进行对比， 相对误差如图7所示。 从图中可以看出， 相对误差随时间的增大而增大， 这是因为随着计算时间的增加， 在t～t+Δt时刻内， 时间t之前累计载荷的松弛效应逐渐增大， 而在计算的过程中忽略了Δt时刻内松弛效应的影响。 虽然相对误差随时间的增大而增大， 但数值仍然较小， 说明了该计算方法可行。 从不同时间增量Δt的计算结果来看， 当Δt从0.001 s缩小10倍后， 相对误差产生了震荡， 并且减小Δt时相对误差降低的幅度不大。 输入条件不变， 拉伸时间改为100 s， Δt取0.001 s， 计算的松弛模量与理论值对比如图8所示。 当拉伸时间为100 s时， 相对误差为0.83%。 采用该方法可以根据复合固体推进剂的拉伸曲线计算松弛模量的变化。
3总结
（1） 由于复合固体推进剂材料的颗粒夹杂效应， 通过仿真计算可知其松弛模量具有明显的应变相关性， 即在不同的定应变水平下测得的试验结果会存在差异， 随着定应变水平的增大， 测得的松弛模量逐渐减小。 计算结果所得结论与züpek[15]试验结果一致。
（2） 根据单胞体模型不同拉伸速率下， 不考虑界面损伤时的拉伸曲线， 计算了拉伸过程中复合固体推进剂松弛模量的变化。 计算结果表明， 在拉伸过程中， 复合固体推进剂的松弛模量有明显的衰减现象， 并且衰减程度与载荷历程相关。
参考文献：
[1] Matous K， Inglis H M， Gu X F， et al. Multiscale Damage Modeling of Solid Propellants： Theory and Computational Framework[C]//41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit， AIAA 2005-4347， 2005： 1-14. [2] Matous K， Geubelle P H. Multiscale Modeling of Particle Debonding in Reinforced Elastomers Subjected to Finite Deformation[J]. International Journal for Numerical Methods in Engineering， 2006， 65（2）： 190-223.
[3] Tan H， Liu C， Huang Y， et al. The Cohesive Law for the Particle/Matrix Interfaces in High Explosives[J]. Journal of the Mechanics and Physics of Solids， 2005， 53（8）： 1892-1917.
[4] Tan H， Huang Y， Liu C， et al. The Uniaxial Tension of Particulate Composite Materials with Nonlinear Interface Debonding[J]. International Journal of Solids and Structures. 2007， 44（6）： 1809-1822.
[5] 史佩， 李高春， 李昊. 复合固体推进剂细观力学模型研究[J]. 计算机仿真， 2007， 24（5）： 21-24.
[6] 赵玖玲， 田先斌. 颗粒复合材料代表性体元并行建模算法研究[J]. 计算机仿真， 2010， 27（1）： 46-49.
[7] 刘著卿， 李高春， 邢耀国. 复合固体推进剂细观损伤扫描电镜实验及数值模拟[J]. 推进技术， 2011， 32（3）： 412-416.
[8] 曲凯， 张旭东， 李高春. 基于内聚力界面脱粘的复合固体推进剂力学性能研究[J]. 火炸药学报， 2008， 31（6）： 77-81.
[9] 李高春， 邢耀国， 王玉峰. 基于细观力学的复合固体推进剂模量预估方法[J]. 推进技术， 2007， 28（4）： 441-444.
[10] 李高春， 邢耀国， 戢治洪， 等. 复合固体推进剂细观界面脱粘有限元分析[J]. 复合材料学报， 2011， 28（3）： 229-235.
[11] 刘承武， 阳建红， 陈飞. 改进的Mori-Tanaka法在复合推进剂非线界面脱粘中的应用[J]. 固体火箭技术， 2011， 34（1）： 67-70.
[12] 马昌兵， 强洪夫， 武文明， 等. 颗粒增强复合材料有效弹性模量预测的多步法[J]. 固体力学学报， 2010， 31（S1）： 12-16.
[13] Zhi Shijun， Sun Bing， Zhang Jianwei. 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.
[14] 职世君， 孙冰， 张建伟. 基于表面粘结损伤的复合固体推进剂细观损伤数值模拟[J]. 推进技术， 2013， 34（2）： 273-279.
[15] Sebnem züpek. Constitutive Equations for Solid Propellants [D]. University of Texas at Austin， 1997.