复合材料学报 第36卷第10期10月2019年 Acta Materiae Compositae sinica Vo.36No.10 Oct 2019 DOI:10.13801/j.cnki. Thclxb.20190221.001 种高颗粒填充率丁羟推进剂二维细观模型 生成方法 张镇国1,侯睫‘1,郜婕1,翁琳2 (1.中国航天科技集团第四研究院第四十一所,西安710025;2.上海工程技术大学城市轨道交通学院,上海201620) 摘要:为了从机制上认识推进剂的宏观力学性能和试验现象,对推进剂开展细观力学研究十分必要,而代表 体积单元(RvE)则是细观力学研究中的重要部分。在随机序列吸附(RSA)方法的基础上,提出将投递区域离散为 背景网格,再进行二次随机投递,可以生成80o%颗粒含量的细观模型,模型生成效率相较于传统RSA方法有所 提高。采用 Python语言编程,利用 ABAQUS的二次开发接口,实现了推进剂RVE模型的快捷生成,并加载周期性 边界条件。采用模型生成实例说明了本文方法的有效性,为开展推进剂的细观分析提供了参考依据 关键词:代表体积单元(RVE)模型;细观力学;背景网格;固体推进剂;周期性边界条件 中图分类号:TB330.1文献标志码:A文章编号:1000-3851(2019)10-2302-06 A method of generating two-dimensional mesoscopic model for hydrox-yl-terminated polybutadiene propellant with high particle volume fraction ZHANG Zhenguo, HOU Xiao, GAO Jie, WENG Lin2 (1. The 41st Institute of the Fourth Academy of CASC, Xi'an 710025, China: 2. School of Urban Railway Transportation, Shanghai University of Engineering Science, Shanghai 201620, China) Abstract: In order to understand the macroscopic mechanical properties and experimental phenomena of propellants from the mechanism, it is necessary to carry out meso-mechanical research on propellants, and the representative volume element(RVE) is an important part of meso-mechanics research. Based on the random sequential adsorption (RSA)method, the delivery area was discretized into a background grid, and then a second random delivery was performed, which can generate a mesoscopic model with 80vol% particle volume fraction, and the model generatio fficiency is also improved compared with the traditional RSA method. The method was programmed in Python lan- guage, and the secondary development interface of ABAQUS was used to realize the rapid generation of RVE model and load periodic boundary conditions. The numerical examples were used to illustrate the effectiveness of the alge rithm, which provides a reference for the meso-analysis of propellants Keywords: representative volume element(RVE) model: meso-mechanical: background grid; solid propellant 目前广泛使用的丁羟( Hydroxyl-terminated粒形状、体积分数、分布规律及各组分相在细观尺 polybutadiene,HTPB)固体推进剂是一种由金属度上的力学响应,对宏观力学行为有着决定性的作 颗粒、氧化剂粒子与具有黏弹性的高分子粘结剂基用-。由于细观力学行为表征固有的复杂性、实验 体复合而成的一种颗粒增强含能复合材料,各组分成本高昂、实验周期长及实验设备、测试方法的局 间有明显的界面,细观结构不均匀,属于非均质材限性,针对细观结构模型的数值仿真技术具有极其 料。复合固体推进剂的细观结构,如颗粒粒度、颗重要的意义[3 收稿日期:2018-10-22;录用日期:2019-01-17;网络出版时间:2019-02-2210:3 网络出版地址:htps://doi.org/10.13801/j. cnki. fhlb.20190221.001 基金项目:固体发动机共性基础研究专项 通讯作者:侯晓,博土,教授,博士生导师,研究方向为药柱结构完整性E-mail:409489678@qg.com 引用格式:张镇国,侯晓,郜婕,等.一种高颗粒填充率丁羟推进剂二维细观模型生成方法[冂].复合材料学报,2019,36(10): ZHANG Z G, HOU X, GAO J, et al. A method of generating two-dimensional mesoscopic model for hydroxyl-termina diene propellant with high particle volume fraction[J]. Acta Materiae Compositae Sinica, 2019, 36(10): 2302-2307(in Chinese)
复合材料学报 第36卷 第10期 10月 2019年 ActaMateriaeCompositaeSinica Vol.36 No.10 Oct 2019 DOI:10.13801/j.cnki.fhclxb.20190221.001 收稿日期:2018-10-22;录用日期:2019-01-17;网络出版时间:2019-02-22 10:32 网络出版地址:https://doi.org/10.13801/j.cnki.fhclxb.20190221.001 基金项目:固体发动机共性基础研究专项 通讯作者:侯晓,博士,教授,博士生导师,研究方向为药柱结构完整性 E-mail:409489678@qq.com 引用格式:张镇国,侯晓,郜婕,等.一种高颗粒填充率丁羟推进剂二维细观模型生成方法[J].复合材料学报,2019,36(10):2302-2307. ZHANGZG,HOUX,GAOJ,etal.Amethodofgeneratingtwo-dimensionalmesoscopicmodelforhydroxyl-terminatedpolybutadienepropellantwithhighparticlevolumefraction[J].ActaMateriaeCompositaeSinica,2019,36(10):2302-2307(inChinese). 一种高颗粒填充率丁羟推进剂二维细观模型 生成方法 张镇国1, 侯晓*1, 郜婕1, 翁琳2 (1.中国航天科技集团第四研究院第四十一所,西安 710025;2.上海工程技术大学 城市轨道交通学院,上海 201620) 摘 要: 为了从机制上认识推进剂的宏观力学性能和试验现象,对推进剂开展细观力学研究十分必要,而代表 体积单元(RVE)则是细观力学研究中的重要部分。在随机序列吸附(RSA)方法的基础上,提出将投递区域离散为 背景网格,再进行二次随机投递,可以生成80vol%颗粒含量的细观模型,模型生成效率相较于传统 RSA 方法有所 提高。采用Python语言编程,利用 ABAQUS的二次开发接口,实现了推进剂 RVE模型的快捷生成,并加载周期性 边界条件。采用模型生成实例说明了本文方法的有效性,为开展推进剂的细观分析提供了参考依据。 关键词: 代表体积单元(RVE)模型;细观力学;背景网格;固体推进剂;周期性边界条件 中图分类号: TB330.1 文献标志码: A 文章编号: 1000-3851(2019)10-2302-06 Amethodofgeneratingtwo-dimensionalmesoscopicmodelforhydrox-yl-terminated polybutadienepropellantwithhighparticlevolumefraction ZHANGZhenguo1,HOU Xiao*1,GAOJie1,WENGLin2 (1.The41stInstituteoftheFourthAcademyofCASC,Xi’an710025,China;2.SchoolofUrbanRailway Transportation,ShanghaiUniversityofEngineeringScience,Shanghai201620,China) Abstract: Inordertounderstandthemacroscopicmechanicalpropertiesandexperimentalphenomenaofpropellants fromthemechanism,itisnecessarytocarryoutmeso-mechanicalresearchonpropellants,andtherepresentative volumeelement(RVE)isanimportantpartofmeso-mechanicsresearch.Basedontherandomsequentialadsorption (RSA)method,thedeliveryareawasdiscretizedintoabackgroundgrid,andthenasecondrandomdeliverywas performed,whichcangenerateamesoscopicmodelwith80vol% particlevolumefraction,andthemodelgeneration efficiencyisalsoimprovedcomparedwiththetraditionalRSAmethod.ThemethodwasprogrammedinPythonlanguage,andthesecondarydevelopmentinterfaceofABAQUSwasusedtorealizetherapidgenerationofRVEmodel andloadperiodicboundaryconditions.Thenumericalexampleswereusedtoillustratetheeffectivenessofthealgorithm,whichprovidesareferenceforthemeso-analysisofpropellants. Keywords: representativevolumeelement(RVE)model;meso-mechanical;backgroundgrid;solidpropellant; periodicboundarycondition 目 前 广 泛 使 用 的 丁 羟 (Hydroxyl-terminated polybutadiene,HTPB)固体推 进 剂 是 一 种 由 金 属 颗粒、氧化剂粒子与具有黏弹性的高分子粘结剂基 体复合而成的一种颗粒增强含能复合材料,各组分 间有明显的界面,细观结构不均匀,属于非均质材 料。复合固体推进剂的细观结构,如颗粒粒度、颗 粒形状、体积分数、分布规律及各组分相在细观尺 度上的力学响应,对宏观力学行为有着决定性的作 用[1-2]。由于细观力学行为表征固有的复杂性、实验 成本高昂、实验周期长及实验设备、测试方法的局 限性,针对细观结构模型的数值仿真技术具有极其 重要的意义[3-5]
张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 细观结构模型的建立是进行细观分析的基础,式中:f为颗粒填充率;V为颗粒体积分数。 复合固体推进剂等复合材料各相的组成虽然在细观1.1颗粒库生成方法 尺度上具有任意性和随机性,但在统计意义上一般 推进剂中的高氯酸铵( Ammonium perchlo 认为具有周期重复性η。因此,人们提出采用一小rate,AP)颗粒总体呈现双峰分布规律,小颗粒 段具有周期性的细观结构来代表整个复合材料,即的平均直径约为20m,大颗粒平均直径在100 代表体积单元( Representative volume element,300pm范围内;A颗粒的粒径范围一般为10 RVE),整个材料看成是RVE周期性排列而成的。30mm。Al颗粒服从正态分布,AP颗粒服从双峰 RVE是关联推进剂细观结构和宏观力学行为的重分布。根据研究需要采用本方法也可以生成只包含 要纽带,具有宏观足够小及细观足够大的特点 AP颗粒的RVE模型。 目前细观结构模型的构建主要基于CT(Com 双峰分布由两个正态分布混合而成。假设X puted tomography)扫描技及几何拓扑方法。N(m,a),Y~N(P2,函),则由X和Y混合而成 马昌兵对HTPB推进剂进行了不受力情况下细的双峰分布为 观结构的微CT扫描,得到了二维断层图像和三维z=LX+(1-1n)Y 空间颗粒分布的重构。CT扫描技术可以真实展现式中,L为权重系数。 推进剂的细观结构,但是重构模型的精度受限于 根据推进剂的级配信息,按照上述正态分布和双 CT扫描设备的精度,且设备价格高昂;几何拓扑峰分布,分别生成各个级配段的颗粒信息(包含颗粒半 方法主要包括随机序列吸附( Random sequential径及颗粒数量),然后按照粒径从大到小的顺序对所 adsorption,RSA)方法[1和分子动力学( Molecu有颗粒进行排序,Al颗粒和AP颗粒分别标记不同 lar dynamics,MD)方法11。RSA方法较简单 的身份信息,组成待投递颗粒库供下一阶段使用。 先生成一定体积分数的颗粒,再在指定的区域内进1.2基于背景网格的颗粒投放方法 行随机投递,但当目标体积分数较高时,传统方法 颗粒从生成到投放完成需要确定3个随机变 需要进行大量的运算,投递效率较低,需要花费大量,即半径r(p)和形心坐标(x(p),y(p)。上述 量的时间;MD方法运算效率较高,可满足高体积生成颗粒库之后,即所有颗粒的半径已经确定,下 含量模型的生成需求,但是MD方法的原理较复面需要进一步确定投放的形心坐标。传统的RSA 杂,编程实现具有很大难度。由于三维的RVE投递方法每次投递时都是在相同的投递区域内生成 模型计算量大且收敛难度高,目前大多学者随机坐标位置,然后判断所投递颗粒与已有颗粒是 主要开展基于二维模型的细观力学研究。 否重叠;随着区域内已投递颗粒的增多,产生有效 为了生成高颗粒填充率(大于65vol%)的HT坐标的难度越来越高。为了提高投放效率,本文将 PB推进剂二维细观RVE模型,本文在传统RSA投递区域首先离散为背景网格,在每一个格点附近 方法的基础上,采用增加背景网格的方法,对投递产生颗粒坐标后便把该格点删除,即将颗粒在区域 区域进行空间离散,减少了无效的投递尝试,且模内服从均匀分布转化为条件均匀分布,进行投递的 型采用周期性几何边界条件。修正后的方法采用具体操作步骤如下所述 首先将投递区域按照一定尺寸离散化,形成如 ABAQUS软件进行后续分析 图1所示的背景网格,将网格点存储为一个 Python 1HIPB推进剂二维细观RvE模型生成方法列表A;然后采用 numpy包中的 random随机函数 目前大多数文献直接采用推进剂真实的颗产生一个0~n之间的随机正整数k,在A中查找 粒填充体积分数作为二维细观模型的目标颗粒填充第k个坐标信息,即对应L4,如下式: 率,赵玖玲等[明通过微CT断面扫描试验证实颗粒A=[L。,L1,L (3) 的二维填充率与其真实颗粒体积分数较吻合。本文L1=(x-1,y-1) (4) 也将采用推进剂的真实颗粒体积含量作为目标填充式中:A表示存储列表;L表示坐标点 率来生成二维细观模型,即 为了进一步增加随机性,实际的颗粒坐标按照 fr (1)下式进一步随机化:
细观结构模型的建立是进行细观分析的基础, 复合固体推进剂等复合材料各相的组成虽然在细观 尺度上具有任意性和随机性,但在统计意义上一般 认为具有周期重复性[6-7]。因此,人们提出采用一小 段具有周期性的细观结构来代表整个复合材料,即 代 表 体 积 单 元 (Representativevolumeelement, RVE),整个材料看成是 RVE 周期性排列而成的。 RVE是关联推进剂细观结构和宏观力学行为的重 要纽带,具有宏观足够小及细观足够大的特点。 目前细观结构模型的构建主要基于 CT(Computedtomography)扫 描 技 术[8]及 几 何 拓 扑 方 法。 马昌兵[9]对 HTPB 推进剂进行了不受力情况下细 观结构的微 CT 扫描,得到了二维断层图像和三维 空间颗粒分布的重构。CT 扫描技术可以真实展现 推进剂的细观结构,但是重构模型的精度受限于 CT 扫描设备的精度,且设备价格高昂;几何拓扑 方法主 要 包 括 随 机 序 列 吸 附 (Random sequential adsorption,RSA)方法[10-13]和分子动力学(Moleculardynamics,MD)方 法[14-16]。RSA 方 法 较 简 单, 先生成一定体积分数的颗粒,再在指定的区域内进 行随机投递,但当目标体积分数较高时,传统方法 需要进行大量的运算,投递效率较低,需要花费大 量的时间;MD方法运算效率较高,可满足高体积 含量模型的生成需求,但是 MD 方法的原理较复 杂[17],编程实现具有很大难度。由于三维的 RVE 模型计算量大且收敛难度高,目前大多学者[3-6,18] 主要开展基于二维模型的细观力学研究。 为了生成高颗粒填充率(大于65vol%)的 HTPB推进剂二维细观 RVE 模型,本文在传统 RSA 方法的基础上,采用增加背景网格的方法,对投递 区域进行空间离散,减少了无效的投递尝试,且模 型采用周期性几何边界条件。修正后的方法采用 Python 语 言 进 行 编 写 实 现, 可 以 方 便 地 导 入 ABAQUS软件进行后续分析。 1 HTPB推进剂二维细观 RVE模型生成方法 目前大多数文献[3-6]直接采用推进剂真实的颗 粒填充体积分数作为二维细观模型的目标颗粒填充 率,赵玖玲等[19]通过微 CT 断面扫描试验证实颗粒 的二维填充率与其真实颗粒体积分数较吻合。本文 也将采用推进剂的真实颗粒体积含量作为目标填充 率来生成二维细观模型,即 fr =Vf (1) 式中:fr 为颗粒填充率;Vf 为颗粒体积分数。 1.1 颗粒库生成方法 推 进 剂 中 的 高 氯 酸 铵 (Ammonium perchlorate,AP)颗粒总体呈现双峰分布规律[20],小颗粒 的平均直径约为20μm,大颗粒平均直径在100~ 300μm 范围内;Al颗粒的粒径范围一般为 10~ 30μm。Al颗粒服从正态分布,AP 颗粒服从双峰 分布。根据研究需要采用本方法也可以生成只包含 AP颗粒的 RVE模型。 双峰分布由两个正态分布混合而成。假设X ~ N(μ1,σ2 1),Y ~ N(μ2,σ2 2),则由 X 和Y 混合而成 的双峰分布为 Z =IaX + (1-Ia)Y (2) 式中,Ia 为权重系数。 根据推进剂的级配信息,按照上述正态分布和双 峰分布,分别生成各个级配段的颗粒信息(包含颗粒半 径及颗粒数量),然后按照粒径从大到小的顺序对所 有颗粒进行排序,Al颗粒和 AP颗粒分别标记不同 的身份信息,组成待投递颗粒库供下一阶段使用。 1.2 基于背景网格的颗粒投放方法 颗粒从生成到投放完成需要确定3个随机变 量,即半径r(p)和形心坐标(x0(p),y0(p))。上述 生成颗粒库之后,即所有颗粒的半径已经确定,下 面需要进一步确定投放的形心坐标。传统的 RSA 投递方法每次投递时都是在相同的投递区域内生成 随机坐标位置,然后判断所投递颗粒与已有颗粒是 否重叠;随着区域内已投递颗粒的增多,产生有效 坐标的难度越来越高。为了提高投放效率,本文将 投递区域首先离散为背景网格,在每一个格点附近 产生颗粒坐标后便把该格点删除,即将颗粒在区域 内服从均匀分布转化为条件均匀分布,进行投递的 具体操作步骤如下所述。 首先将投递区域按照一定尺寸离散化,形成如 图1所示的背景网格,将网格点存储为一个Python 列表A;然后采用numpy包中的random 随机函数 产生一个0~n之间的随机正整数k,在 A 中查找 第k 个坐标信息,即对应Lk-1,如下式: A = [L0,L1,L2,…,Ln] (3) Lk-1 = (xk-1,yk-1) (4) 式中:A 表示存储列表;L 表示坐标点。 为了进一步增加随机性,实际的颗粒坐标按照 下式进一步随机化: 张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 ·2303·
2304 复合材料学报 h—2h +hR 特征尺寸的附近;当模型中颗粒填充率低于68vol% 5)时,可以适当选取较大的网格尺寸。背景网格尺寸越 小,则程序运行时间越长,但生成模型的相对误差会 越小 式中,R为 numpy包中 random函数生成的0~ 之间均匀分布的随机数 确定实际的投递坐标点之后,通过下式判断其 图1为背景网格投放原理示意图。假设L4-1 是否与已投递的颗粒重叠 为0点,则实际颗粒的位置按照式(5)计算,其可能(1-x)2+(y-y)≥A(m+n) 位于图1中的abcd区域内;需要注意的是当颗粒 圆心到投递区域边界的距离小于颗粒半径时,则需式中,表示影响系数,一般取1.02-1.05。 要在对边对应位置投递相同的颗粒,以保证几何周 影响系数控制着颗粒间的最小间距,当颗粒较 期性边界条件,下式说明了这种变换关系: 多时应适当取小值。当按照上述步骤成功投递完颗 xk-1=xk-1±Hx (6) 粒之后,把该颗粒代表的背景网格信息即Lk1从存储 y1=y1±H (7) 列表中删除,从而在投递下一次颗粒时,其对应的投 式中:x1、y1.为对应位置的新坐标;x4-1、y-1为递区域就不包含该背景网格,从而提高了下次投递 原位置坐标;H,为投递区域x方向长度;H,为投的效率。按照上述步骤依次投递颗粒库的每个颗粒 递区域y方向长度;当原位置坐标靠近投递初始边直至全部投递完成 界xm或ymn时,式(6)和式(7)选择正号,当原位2HIPB推进剂二维细观RⅤE模型生成及数 置坐标靠近投递外缘边界xm或ym时,式(6)和式值仿真实例 (7)选择负号,效果如图1所示。 2.1HIPB推进剂二维细观RⅤE模型生成实例 为了确定背景网格尺寸的合适大小,特别地, 文献[21]中的典型HTPB推进剂配方如表1所 我们假设一个颗粒其圆心正好位于图中c点,且其示,该配方中AP的体积分数为63.8%,其中AP 半径为颗粒库最小半径,另外一个颗粒恰好落在a 大颗粒约占33.8vol%,AP小颗粒约占30vol%,A 点,显然当两个圆恰好相切时,对应的网格尺寸正颗粒占12.35%,颗粒总体积分数约为76.1% 好不遗漏任何可以投递的区域,此时的背景网格尺 寸称为特征网格尺寸,其值满足下式: 表1丁羟(HIPB推进剂的典型配方 lable 1 Typical formulation of hydroxyl-terminated h=√2R (8) polybutadiene(HTPB) propellant 式中:h代表背景网格尺寸;Rm为颗粒库的最小颗 Component HTPB 粒半径。 Mass fraction/wt% 8 69.5 18.5 充率大于68v01%时,网格尺寸的合适大小应在此Noe: AP- Ammonium perchlorate 经过大量的数值模拟试验,发现当模型中颗粒 12.3 Density/(g·cm-3)0.9 填 文献[20]中给出了二级配推进剂的粒径范围 本文按照表2所示的粒径分布参数生成细观模型。 RVE模型的边长选为1mm×1mm;背景网格尺寸h 选为7pm;影响系数λ选为1.05。在 ABAQUS中生 成的考虑小颗粒的颗粒填充HTPB推进剂RVE模 型的实际颗粒填充率如图2所示,实际模型的颗粒 填充率为75.7vol%,相对于目标填充率76.1vol% 误差仅为0.53%,该误差可通过调节影响系数λ进 一步修正 可以看到,如果考虑细小颗粒建立全级配模型 图1背景网格投放原理示意图 则要求必须具有较高的计算资源才可以进行后续分 wagram 析。为了减少网格数量进而节省计算资源,同时也为
xt k-1 =xk-1 - h 2 +hR y t k-1 =yk-1 - h 2 +hR (5) 式中,R 为 numpy包中random 函数生成的0~1 之间均匀分布的随机数。 图1 背景网格投放原理示意图 Fig.1 Schematicdiagramofbackgroundgridplacement 图1为背景网格投放原理示意图。假设 Lk-1 为o点,则实际颗粒的位置按照式(5)计算,其可能 位于图1中的abcd 区域内;需要注意的是当颗粒 圆心到投递区域边界的距离小于颗粒半径时,则需 要在对边对应位置投递相同的颗粒,以保证几何周 期性边界条件,下式说明了这种变换关系: x'k-1=xt k-1 ±Hx (6) y'k-1=y t k-1 ±Hy (7) 式中:x'k-1、y'k-1 为对应位置的新坐标;xt k-1、y t k-1 为 原位置坐标;Hx 为投递区域x 方向长度;Hy 为投 递区域y 方向长度;当原位置坐标靠近投递初始边 界xmin或ymin时,式(6)和式(7)选择正号,当原位 置坐标靠近投递外缘边界xmax或ymax时,式(6)和式 (7)选择负号,效果如图1所示。 为了确定背景网格尺寸的合适大小,特别地, 我们假设一个颗粒其圆心正好位于图中c点,且其 半径为颗粒库最小半径,另外一个颗粒恰好落在a 点,显然当两个圆恰好相切时,对应的网格尺寸正 好不遗漏任何可以投递的区域,此时的背景网格尺 寸称为特征网格尺寸,其值满足下式: h= 2Rmin (8) 式中:h代表背景网格尺寸;Rmin为颗粒库的最小颗 粒半径。 经过大量的数值模拟试验,发现当模型中颗粒 填充率大于68vol%时,网格尺寸的合适大小应在此 特征尺寸的附近;当模型中颗粒填充率低于68vol% 时,可以适当选取较大的网格尺寸。背景网格尺寸越 小,则程序运行时间越长,但生成模型的相对误差会 越小。 确定实际的投递坐标点之后,通过下式判断其 是否与已投递的颗粒重叠: (xt k-1 -xi)2 +(y t k-1 -y i)2 ≥λ(rt +ri) i=0,1,2,…,t-1 (9) 式中,λ表示影响系数,一般取1.02~1.05。 影响系数控制着颗粒间的最小间距,当颗粒较 多时应适当取小值。当按照上述步骤成功投递完颗 粒之后,把该颗粒代表的背景网格信息即Lk-1从存储 列表中删除,从而在投递下一次颗粒时,其对应的投 递区域就不包含该背景网格,从而提高了下次投递 的效率。按照上述步骤依次投递颗粒库的每个颗粒 直至全部投递完成。 2 HTPB推进剂二维细观 RVE模型生成及数 值仿真实例 2.1 HTPB推进剂二维细观 RVE模型生成实例 文献[21]中的典型 HTPB推进剂配方如表1所 示,该配方中 AP的体积分数为63.8vol%,其中 AP 大颗粒约占33.8vol%,AP小颗粒约占30vol%,Al 颗粒占12.35%,颗粒总体积分数约为76.1%。 表1 丁羟(HTPB)推进剂的典型配方[21] Table1 Typicalformulationofhydroxyl-terminated polybutadiene(HTPB)propellant [21] Component HTPB AP Al Other Massfraction/wt% 8.0 69.5 18.5 4 Volume fraction/vol% 23.9 63.8 12.3 — Density/(g·cm-3) 0.9 1.95 2.7 — Notes:AP—Ammoniumperchlorate. 文献[20]中给出了二级配推进剂的粒径范围, 本文按照表2所示的粒径分布参数生成细观模型。 RVE模型的边长选为1mm×1mm;背景网格尺寸h 选为7μm;影响系数λ选为1.05。在 ABAQUS中生 成的考虑小颗粒的颗粒填充 HTPB推进剂 RVE模 型的实际颗粒填充率如图2所示,实际模型的颗粒 填充率为75.7vol%,相对于目标填充率76.1vol%, 误差仅为0.53%,该误差可通过调节影响系数λ进 一步修正。 可以看到,如果考虑细小颗粒建立全级配模型, 则要求必须具有较高的计算资源才可以进行后续分 析。为了减少网格数量进而节省计算资源,同时也为 ·2304· 复合材料学报
张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 了降低收敛难度,很多文献先根据细观均匀化理 论如Mori- Tanaka方法将小粒径的Al颗粒、AP颗 粒与HTPB基体进行均匀化操作,得到均匀化之后 的新参数作为基体参数,细观模型则只考虑大粒径 的AP。只考虑表1推进剂配方中的粗AP颗粒,生 成的细观RVE模型如图3所示,实际模型的填充率 为34.5vol%,相对于目标填充率33.8vo1%,误差为 2.1%;图3(b)所示图形为将所生成模型阵列之后的 a)Rve model generated (b)Effect of arraying by ABAQUS 效果,可见其几何边界满足周期性边界条件。 图3只考虑大颗粒的颗粒填充HTPB推进剂RVE模型的 表2二级配HTPB推进剂的粒径分布参数 实际颗粒填充率(34.7vol%) Table 2 Particle size distribution parameters of Fig 3 Actual filling fraction (34 7vol %)of RVE model of tworgraded HTPB propellant particle filling hTPB propellant considering only large particles Component Distribution parameter-(.02) 表3五种配方的HTPB推进剂各组分的体积含量 Thick ap Table 3 Volume fraction of each component of five HTPB propellants Thin Ap Formula HTPB/vol% Thick AP/vol% Thin AP/vol% Al/vol% Notes: /eXpectation value: d2-Variance: d-Standard deviation. C 40 20 图4为本文算法与传统RSA算法计算的HT PB推进剂RVE模型生成效率比较。可以看到,在 颗粒总体积分数小于65vol%时,两种方法生成模 型所需要的时间比较接近;当颗粒总体积分数大于 65vol%时,采用传统的RSA方法所需要的运行时 间就开始快速增加,采用本文的背景网格方法所需 要的运行时间增加则相对缓慢,对于生成80vol% 颗粒体积含量的细观模型,本文的修正方法运行时 图2考虑小颗粒的颗粒填充HTPB推进剂RVE模型的 实际颗粒填充率(75.7vol%) Fig 2 Actual filling fraction(75. 7 vol %)of RVE model 1 000F+-Relative error of traditional RSA algorithms/A article filling HTPB propellant considering small particles -Relative error of algorithms in this papey 0.03 2.2HTPB推进剂二维细观RⅤE模型生成效率 为了对比采用背景网格方法与传统RSA方法 0.02 的效率,本文分别采用背景网格方法和传统的RSA220 0.01 方法来生成五种配方的推进剂RVE模型,各组分 的具体体积分数见表3。配方A的颗粒总体积含量 最小为60vol%,配方E的颗粒总体积含量最高为 55 70 58085 80vo1%;RVE的大小均设置为1mm×1mm,控 图4本文算法与传统RSA算法HTPB推进剂RVE模型 制生成模型的实际填充率的误差在5%以内,并且 生成效率的比较 为了减小随机性的影响,每个配方分别采用两种方 Fig 4 Comparison of the generation efficiency of HTPB 法各生成十次符合误差要求的模型,最后对该配方 propellant RvE model by the algorithms in this paper 各个方法的十次结果求算数平均作为对比依据 with traditional RSa algorithm
了降低收敛难度,很多文献[3-6]先根据细观均匀化理 论如 Mori-Tanaka 方法将小粒径的 Al颗粒、AP颗 粒与 HTPB基体进行均匀化操作,得到均匀化之后 的新参数作为基体参数,细观模型则只考虑大粒径 的 AP。只考虑表1推进剂配方中的粗 AP颗粒,生 成的细观 RVE模型如图3所示,实际模型的填充率 为34.5vol%,相对于目标填充率33.8vol%,误差为 2.1%;图3(b)所示图形为将所生成模型阵列之后的 效果,可见其几何边界满足周期性边界条件。 表2 二级配 HTPB推进剂的粒径分布参数 Table2 Particlesizedistributionparametersof two-gradedHTPBpropellant Component Distributionparameterr~ N(μ,σ2) ThickAP μ=90μm σ=25μm ThinAP μ=10μm σ=8μm Al μ=10μm σ=8μm Notes:μ—Expectationvalue;σ2—Variance;σ—Standarddeviation. 图2 考虑小颗粒的颗粒填充 HTPB推进剂 RVE模型的 实际颗粒填充率(75.7vol%) Fig.2 Actualfillingfraction(75.7vol%)ofRVEmodelof particlefillingHTPBpropellantconsideringsmallparticles 2.2 HTPB推进剂二维细观 RVE模型生成效率 为了对比采用背景网格方法与传统 RSA 方法 的效率,本文分别采用背景网格方法和传统的RSA 方法来生成五种配方的推进剂 RVE 模型,各组分 的具体体积分数见表3。配方 A 的颗粒总体积含量 最小为60vol%,配方 E 的颗粒总体积含量最高为 80vol%;RVE的大小均设置为1mm×1mm,控 制生成模型的实际填充率的误差在5%以内,并且 为了减小随机性的影响,每个配方分别采用两种方 法各生成十次符合误差要求的模型,最后对该配方 各个方法的十次结果求算数平均作为对比依据。 图3 只考虑大颗粒的颗粒填充 HTPB推进剂 RVE模型的 实际颗粒填充率(34.7vol%) Fig.3 Actualfillingfraction(34.7vol%)ofRVEmodelof particlefillingHTPBpropellantconsideringonlylargeparticles 表3 五种配方的 HTPB推进剂各组分的体积含量 Table3 Volumefractionofeachcomponentof fiveHTPBpropellants Formula HTPB/vol% ThickAP/vol% ThinAP/vol% Al/vol% A 40 35 15 10 B 35 35 20 10 C 30 40 20 10 D 24 40 24 12 E 20 40 25 15 图4 本文算法与传统 RSA 算法 HTPB推进剂 RVE模型 生成效率的比较 Fig.4 ComparisonofthegenerationefficiencyofHTPB propellantRVEmodelbythealgorithmsinthispaper withtraditionalRSAalgorithm 图4为本文算法与传统 RSA 算法计算的 HTPB推进剂 RVE模型生成效率比较。可以看到,在 颗粒总体积分数小于65vol%时,两种方法生成模 型所需要的时间比较接近;当颗粒总体积分数大于 65vol%时,采用传统的 RSA 方法所需要的运行时 间就开始快速增加,采用本文的背景网格方法所需 要的运行时间增加则相对缓慢,对于生成80vol% 颗粒体积含量的细观模型,本文的修正方法运行时 张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 ·2305·
复合材料学报 间大约是传统RSA运行时间的一半;这是由于在 Stress/MPa 目标填充率较高时,传统RSA方法运行后期需要(Awg:75%) +7,565e+00 进行大量的投递尝试才可以成功生成一个有效的投 递位置,而增加背景网格后,当该格点附近投递完 +3900 +8.000e-01 颗粒后,该格点即删除,一定程度上减少了无效的 投递尝试,使生成效率有所提升;并且可以看到采 用传统RSA方法生成模型的相对误差整体大于背 景网格方法的相对误差。综合分析可知本文提出的 增加背景网格的方法在生成高颗粒填充率(大于 65vol%)的细观模型时相较传统RSA方法具有 定的优势。 2.3HIPB推进剂二维细观RVE模型仿真实例 考虑到计算资源有限,本文选用2.1节生成的 只考虑粗AP的简化颗粒模型进行仿真验证。通过 图5应变e=0.15时HTPB推进剂 von mises应力分布 Fig 5 Von Mises stress distribution of HTPB propellant ABAQUS的多点约束( Multi-point constraints, at strain of 0. 15 MPC)功能定义周期性边界条件[21,Y向施加 0.3mm的位移载荷;采用文献[3]中的材料参数 孔洞,小颗粒处的脱湿现象仍不太明显;已经脱湿 AP颗粒采用线弹性本构;本文只考虑复合基体的 的颗粒之间的基体部分 Mises应力较大,预示着空 超弹特性,忽略其黏弹特性,复合基体采用Ne 穴之间容易发生贯穿,从而形成明显的裂纹导致断 Hookean超弹性本构,根据文献[3]的试验数据拟 裂失效。 合得到本构参数;采用双线性内聚力模型模拟界面3结论 分离,目前还没有有效的获取界面参数的方法,本 (1)针对丁羟(HTPB)推进剂,本文通过将颗 文参考文献[6]及根据经验选取界面参数;HTPB粒投放区域离散化为背景网格,提高了颗粒投放效 推进剂的弹性模量(E)、泊松比(ν)及基体的Neo-率,弥补了传统的随机序列吸附(RSA)方法在生成 Hookean本构参数Cuo如表4所示,HTPB/AP界高颗粒含量(大于65vol%)模型时效率低下的不足 面参数如表5所示 且方法可以拓展到三维的情况 表4简化配方中HTPB推进剂中各组分的力学属性 (2)数值分析算例表明本文所建立的HTPB Table 4 Mechanical properties of components in 推进剂代表体积单元(RVE)的有效性,推进剂的细 simplified formulation of HTPB propellant 观典型脱湿特征可以完整地模拟出来,为推进剂的 E/MPa Co/MPa 细观模拟提供了一定参考。 32447 0.1433 (3)方法采用 Python语言进行编写,可以方 0.4991 Notes: E- elasticity modul= Poissovis ratio: cia-c0msm便地利用 ABAQUS的二次开发接口将生成的几何 模型导入 ABAQUS进行后续分析,为推进剂的细 表5HTPB/AP界面损伤参数 观RVE模型的建立提供了一种简单高效的方式 Table 5 Damage parameters of HTPB/AP interface 参考文献 Initial damag racture energy/ stress/MPa (MPa.mm-) 1]陈雄,许进升,郑健.固体推进剂黏弹性力学[M].北京:北 1000 京理工大学出版社,2016:21-22 数值计算结果如图5所示。可以看到,RVE模 HEN X, XU JS, ZHENG J. Viscoelastic mechanics of sol- 型的边界处满足位移连续和应力连续条件,且可以 id propellantsLMI. Beijing: Beijing Institute of Technology Press, 2016:21-22 (in Chinese) 模拟出典型的脱湿特征,并且从计算结果来看,脱2]刘新国,刘佩进,强洪夫,复合固体推进剂脱湿研究进 湿首先发生在大颗粒处,如图5所示,在应变为 LJ].固体火箭技术,2018,41(3):313-318 0.15时,大颗粒处界面已经明显失效并形成了微 LIU X G, LIU PJ, QIANG H F. Recent advances
间大约是传统 RSA 运行时间的一半;这是由于在 目标填充率较高时,传统 RSA 方法运行后期需要 进行大量的投递尝试才可以成功生成一个有效的投 递位置,而增加背景网格后,当该格点附近投递完 颗粒后,该格点即删除,一定程度上减少了无效的 投递尝试,使生成效率有所提升;并且可以看到采 用传统 RSA 方法生成模型的相对误差整体大于背 景网格方法的相对误差。综合分析可知本文提出的 增加背 景 网 格 的 方 法 在 生 成 高 颗 粒 填 充 率 (大 于 65vol%)的细观模型时相较传统 RSA 方法具有一 定的优势。 2.3 HTPB推进剂二维细观 RVE模型仿真实例 考虑到计算资源有限,本文选用2.1节生成的 只考虑粗 AP的简化颗粒模型进行仿真验证。通过 ABAQUS 的 多 点 约 束 (Multi-pointconstraints, MPC)功 能 定 义 周 期 性 边 界 条 件[22],Y 向 施 加 0.3mm 的位移载荷;采用文献[3]中的材料参数, AP颗粒采用线弹性本构;本文只考虑复合基体的 超弹特性,忽略其黏弹特性,复合基体采用 NeoHookean超弹性本构,根据文献[3]的试验数据拟 合得到本构参数;采用双线性内聚力模型模拟界面 分离,目前还没有有效的获取界面参数的方法,本 文参考文献[6]及根据经验选取界面参数;HTPB 推进剂的弹性模量(E)、泊松比(ν)及基体的 NeoHookean本构参数 C10 如表4所示,HTPB/AP 界 面参数如表5所示。 表4 简化配方中 HTPB推进剂中各组分的力学属性[3] Table4 Mechanicalpropertiesofcomponentsin simplifiedformulationofHTPBpropellant Component E/MPa ν C10/MPa AP 32447 0.1433 — Compositematrix 7.393 0.4991 1.2328 Notes:E—Elasticity modulus;ν—Poissovisratio;C10—Constitutiveparameter. 表5 HTPB/AP界面损伤参数 Table5 DamageparametersofHTPB/APinterface Initialdamage stress/MPa Initialstiffness/ (MPa·mm-1) Fractureenergy/ (J·m-2) 0.5 1000 6 数值计算结果如图5所示。可以看到,RVE模 型的边界处满足位移连续和应力连续条件,且可以 模拟出典型的脱湿特征,并且从计算结果来看,脱 湿首先发生在大颗粒处,如图 5 所示,在应变为 0.15时,大颗粒处界面已经明显失效并形成了微 图5 应变ε=0.15时 HTPB推进剂vonMises应力分布 Fig.5 VonMisesstressdistributionofHTPBpropellant atstrainof0.15 孔洞,小颗粒处的脱湿现象仍不太明显;已经脱湿 的颗粒之间的基体部分 Mises应力较大,预示着空 穴之间容易发生贯穿,从而形成明显的裂纹导致断 裂失效。 3 结 论 (1)针对丁羟(HTPB)推进剂,本文通过将颗 粒投放区域离散化为背景网格,提高了颗粒投放效 率,弥补了传统的随机序列吸附(RSA)方法在生成 高颗粒含量(大于65vol%)模型时效率低下的不足, 且方法可以拓展到三维的情况。 (2)数值分析算例表明本文所建立的 HTPB 推进剂代表体积单元(RVE)的有效性,推进剂的细 观典型脱湿特征可以完整地模拟出来,为推进剂的 细观模拟提供了一定参考。 (3)方法采用 Python语言进行编写,可以方 便地利用 ABAQUS的二次开发接口将生成的几何 模型导入 ABAQUS进行后续分析,为推进剂的细 观 RVE模型的建立提供了一种简单高效的方式。 参考文献: [1] 陈雄,许进升,郑健.固体推进剂黏弹性力学[M].北京:北 京理工大学出版社,2016:21-22. CHENX,XUJS,ZHENGJ.Viscoelasticmechanicsofsolidpropellants[M].Beijing:BeijingInstituteofTechnology Press,2016:21-22(inChinese). [2] 刘新国,刘 佩 进,强 洪 夫.复 合 固 体 推 进 剂 脱 湿 研 究 进 展 [J].固体火箭技术,2018,41(3):313-318. LIU X G,LIU P J,QIANG H F.Recentadvanceson ·2306· 复合材料学报
张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 2307 research of the dewetting in composite solid propellants [J]. unit cell models for metal randomly oriented discontinuou Journal of Solid Rocket Technology, 2018, 41(3): 313-318 reinforcementsLJI. Computational Matrix Composites Mate- (in Chinese) rials science,2005,25(1-2):42-53. 3]常武军.复合固体推进剂细观损伤及其数值仿真研究[D].[13]舒刚.超高性能混凝土钢纤维诱蚀的细观力学分析方法研究 南京:南京理工大学,2013 D].成都:西南交通大学,2017 CHANG WJ. Research on microstructural damage and its SHU G. Study on micromechanical analysis method of ultra- numerical simulation method forcomposite solid propellant high performance concrete steel fiber corrosion [D]. Cheng- LDI. Nanjing: Nanjing University of Science and Technology du Southwest Jiaotong University, 2017(in Chinese). 2013(in Chinese) [14 KNOTT G M, JACKSON T L, BUCKMASTER J. Random [4 MATOUS K, INGLIS H M, GUX, et al. Multiscale model packing of heterogeneous propellants [J]. AlAA Journal ing of solid propellants: From particle packing to failure[J]. 2001,39(4):678-686. Composite Science and Technology, 2007, 67 (7-8) [15] KOCHEVETS S V. Random sphere packing model of hetero- geneous propellantsLDI. Champaign: University of Illinois at [5 INGLIS H M, GEUBELLE P H, MATOUS K, et al. Cohe- Urbana Champaign, 2002. ive modeling of dewetting in particle composites Microme- [16] DONEv A, TORQUATO S, STILLINGER F H. Neighbor multiscale finite element analysis [J]. Materials list collision-driven molecular d and Structures, 2007, 39(6):580-595 spherical hard particles[MI. San Diego: Academic Press P 6]职世军,曹付齐,申志彬,等.复合固体推进剂颗粒脱湿损 fessional, Inc, 2005 伤参数反演[J].推进技术,2016,37(10):1977-1983 [17 HAN K, FENG YT, OWEN DR J. Sphere packing with a ZHI SJ, CAO F Q, SHEN Z B, et al. Parameters inversion geometric based compression algorithm[J]. Powder Technol- of particle dewetting damage of composite solid propellan gy,2005,155:33-41 ]. Journal of propulsion Technology,2016,37(10):1977[18]胡大琳,张立兴,陈定市.二维细观随机混凝土模型的建立 983(in Chinese) 和应用[刀].长安大学学报,2017,37(3):53-63 [7]刘著卿,李高春,邢耀国,等.复合固体推进剂细观损伤扫 HU DL, ZHANG L X, CHEN D S. Establishment and 描电镜实验及数值模拟[J].推进技术,2011,32(3) application of 2D mesoscopic stochastic concrete modelLJ] 412-4 Journal of Chang'an University, 2017, 37(3): 53-63(in Chi- LIU Z Q, LI G C, XING Y G, et al. Numerical simulation and SEM study on the microstructural damage of composite[19]赵玖玲,强洪夫.复合固体推进剂宏细观损伤机理[M].北 olid propelantsLJ. Journal of Propulsion Technology 京:中国宇航出版社,2014:114-115 2011,32(3):412-416( in Chinese) ZHAO J L, QIANG H F. Macroscopic and microscopic dam- [8] HU J, QIAN Z D, LIU Y, et al. High-temperature failure in age mechanism of composite solid propellant [M]. Beijing asphalt mixtures using micro-structural investigation and China Aerospace Publishing House, 2014: 114-115 (in Chi- mage analysis [J]. Construction and Building Materials 2015,84:136-145. [20 MATOUS K, GEUBELLE P H. Finite element formulation 9]马昌兵.复合固体推进剂细观结构建模极其力学行为数值模 for modeling particle debonding unreinforced elastomers sub- 拟研究[D]。西安:第二炮兵工程学院,2011 jected to finite deformations [J]. Computer Methods in MA CB Research sostructure modeling and mechani- Applied Mechanics and Engineering, 2006, 196(1-3): cal behaviors numerical simulation of composite solid propel- 620-366. lantLD].Xi’an: The Second Artillery Engineering Institute,[21]侯林法.复合固体推进剂[M].北京:中国宇航出版社 01 (in Chinese) 1994:10-11 10 SEGURADO J, LLORCA J. A new three-dimensional inte HOU L F. Composite solid propellant[M]. Beijing: China face finite element to simulate fracture in compositesLJJ. In- Astronautic Publishing House, 1994: 10-11 (in Chinese) al Journal of Solids&. Structures,2004,41(11):[22]张超,许希武,严雪.纺织复合材料细观力学分析的一般性 2977-2993. 周期性边界条件及其有限元实现[J].航空学报,2012,3 [11] SEGURADO J, LLORCA J. A computational micromechan- (7):1636-1645 ics study of the effect of interface decohesive on the mechani- ZHANG C, XU W, YAN X. General periodic boundary cal behavior of composites[J]. Acta Materialia, 2005 conditions and its application on micromechanical finite ele 4931-4942. ment analysis of textile composites[J]. Acta Aeronautica et [12] BOHM H J, ECKSCHLAGER A, HAN W. Multi-inclusion Astronautica sinica, 2012, 33(7):1636-1645(in Chinese)
researchofthedewettingincompositesolidpropellants[J]. JournalofSolidRocketTechnology,2018,41(3):313-318 (inChinese). [3] 常武军.复合固体推进剂细观损伤及其数值仿真研究[D]. 南京:南京理工大学,2013. CHANG W J.Researchon microstructuraldamageandits numericalsimulation methodforcompositesolid propellant [D].Nanjing:NanjingUniversityofScienceandTechnology, 2013(inChinese). [4] MATOUSK,INGLISH M,GUX,etal.Multiscalemodelingofsolidpropellants:Fromparticlepackingtofailure[J]. Composite Science and Technology, 2007, 67 (7-8): 1694-1708. [5] INGLISH M,GEUBELLEPH,MATOUSK,etal.Cohesivemodelingofdewettinginparticlecomposites Micromechanicsvs.multiscalefiniteelementanalysis[J].Materials andStructures,2007,39(6):580-595. [6] 职世军,曹付齐,申志彬,等.复合固体推进剂颗粒脱湿损 伤参数反演[J].推进技术,2016,37(10):1977-1983. ZHISJ,CAOFQ,SHENZB,etal.Parametersinversion ofparticledewettingdamageofcompositesolidpropellants [J].JournalofPropulsionTechnology,2016,37(10):1977- 1983(inChinese). [7] 刘著卿,李高春,邢耀国,等.复合固体推进剂细观损伤扫 描电 镜 实 验 及 数 值 模 拟 [J].推 进 技 术,2011,32(3): 412-416. LIUZQ,LIGC,XING Y G,etal.Numericalsimulation andSEMstudyonthemicrostructuraldamageofcomposite solid propelants[J].Journal of Propulsion Technology, 2011,32(3):412-416(inChinese). [8] HUJ,QIANZD,LIUY,etal.High-temperaturefailurein asphalt mixtures using micro-structuralinvestigation and imageanalysis[J].Construction and Building Materials, 2015,84:136-145. [9] 马昌兵.复合固体推进剂细观结构建模极其力学行为数值模 拟研究[D].西安:第二炮兵工程学院,2011. MACB.Researchonmesostructuremodelingandmechanicalbehaviorsnumericalsimulationofcompositesolidpropellant[D].Xi’an:TheSecondArtilleryEngineeringInstitute, 2011(inChinese). [10] SEGURADOJ,LLORCAJ.Anewthree-dimensionalinterfacefiniteelementtosimulatefractureincomposites[J].InternationalJournalofSolids& Structures,2004,41(11): 2977-2993. [11] SEGURADOJ,LLORCAJ.Acomputationalmicromechanicsstudyoftheeffectofinterfacedecohesiveonthemechanicalbehaviorofcomposites[J].Acta Materialia,2005,53: 4931-4942. [12] BOHM HJ,ECKSCHLAGERA,HAN W.Multi-inclusion unitcellmodelsfor metalrandomlyorienteddiscontinuous reinforcements[J].ComputationalMatrixComposites MaterialsScience,2005,25(1-2):42-53. [13] 舒刚.超高性能混凝土钢纤维诱蚀的细观力学分析方法研究 [D].成都:西南交通大学,2017. SHU G.Studyonmicromechanicalanalysismethodofultrahighperformanceconcretesteelfibercorrosion[D].Chengdu:SouthwestJiaotongUniversity,2017(inChinese). [14] KNOTTG M,JACKSONTL,BUCKMASTERJ.Random packingofheterogeneouspropellants[J].AIAA Journal, 2001,39(4):678-686. [15] KOCHEVETSSV.Randomspherepackingmodelofheterogeneouspropellants[D].Champaign:UniversityofIllinoisat UrbanaChampaign,2002. [16] DONEV A,TORQUATOS,STILLINGERFH.Neighbor listcollision-driven moleculardynamicssimulationfornonsphericalhardparticles[M].SanDiego:AcademicPressProfessional,Inc,2005. [17] HAN K,FENGYT,OWENDRJ.Spherepackingwitha geometricbasedcompressionalgorithm[J].PowderTechnology,2005,155:33-41. [18] 胡大琳,张立兴,陈定市.二维细观随机混凝土模型的建立 和应用[J].长安大学学报,2017,37(3):53-63. HU D L,ZHANG L X,CHEN D S.Establishmentand applicationof2D mesoscopicstochasticconcrete model[J]. JournalofChang’anUniversity,2017,37(3):53-63(inChinese). [19] 赵玖玲,强洪夫.复合固体推进剂宏细观损伤机理[M].北 京:中国宇航出版社,2014:114-115. ZHAOJL,QIANG HF.Macroscopicandmicroscopicdamagemechanism ofcompositesolidpropellant[M].Beijing: ChinaAerospacePublishing House,2014:114-115 (inChinese). [20] MATOUSK,GEUBELLEP H.Finiteelementformulation formodelingparticledebondinginreinforcedelastomerssubjectedto finite deformations[J].Computer Methodsin Applied Mechanics and Engineering,2006,196 (1-3): 620-366. [21] 侯 林 法.复 合 固 体 推 进 剂 [M].北 京:中 国 宇 航 出 版 社, 1994:10-11. HOULF.Compositesolidpropellant[M].Beijing:China AstronauticPublishingHouse,1994:10-11(inChinese). [22] 张超,许希武,严雪.纺织复合材料细观力学分析的一般性 周期性边界条件及其有限元实现[J].航空学报,2012,33 (7):1636-1645. ZHANGC,XU X W,YAN X.Generalperiodicboundary conditionsanditsapplicationon micromechanicalfiniteelementanalysisoftextilecomposites[J].ActaAeronauticaet AstronauticaSinica,2012,33(7):1636-1645(inChinese). 张镇国,等:一种高颗粒填充率丁羟推进剂二维细观模型生成方法 ·2307·